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FOREWORD 

The  work  reported  here  is  the  result  of  a  continuing  collaboration  between  the  Centre  National 
D'Etudes  Des  Telecommunications,  France,  and  the  National  Bureau  of  Standards. 

From  October  1978  through  September  1979,  Dr.  Michele  Guillaume  was  a  Guest  Scientist  in  the 
Time  Domain  Metrology  Group  of  the  Electromagnetic  Technology  Division,  NBS ,  Boulder  Laboratories,  Boulder, 
Colorado.   At  that  time  Drs.  Guillaume  and  Nahman  worked  together  to  produce  the  experimental  results 
reported  here.   It  has  taken  some  time  for  the  authors  to  complete  their  analyses  of  the  experimental 
data,  to  provide  some  analytic  bases  for  understanding  the  deconvolution  and  filtering  processes,  and  to 
prepare  this  manuscript.   All  of  these  latter  activities  were  accomplished  with  the  authors  located  in 
France  and  the  USA,  respectively.   Thus,  there  has  been  some  delay  in  publishing  this  work. 


DECONVOLUTION  OF  TIME  DOMAIN  WAVEFORMS  IN  THE  PRESENCE  OF  NOISE 

N.  S.  Nahman  and  M.  E.  Guillaume  t 
Electromagnetic  Technology  Division 
National  Bureau  of  Standards 
Boulder,  Colorado  80303 


Deconvolution  or  inverse  filtering  was  used  to  determine  the  impulse  response 
of  a  system  using  noisy  input  and  output  time  domain  sequences  (discrete  data) . 
Frequency  and  time  domain  methods  were  studied  along  with  the  synthesis  of  the 
filters  required  to  obtain  stable  and  smooth  results.   For  the  methods  studied  it 
was  concluded  that  the  superior  technique  was  provided  by  an  optimal  frequency 
domain  method  implemented  via  the  FFT.   Also,  it  is  pointed  out  that  the  time 
domain  methods  are  only  in  their  infancy  and  still  retain  the  promise  of  avoiding 
transform  domain  filtering.   Examples  are  presented  in  which  the  impulse  responses 
are  determined  in  the  presence  of  varying  degrees  of  noise  for  a  coaxial  transmission 
line,  a  wave-shaping  filter,  and  a  broad-band  antenna. 

Key  words:   deconvolution;  impulse  response^  inverse  filtering;  noise  contamination; 
time  domain;  waveforms. 


1.   INTRODUCTION 

The  analogy  between  a  physical  linear  system  and  the  convolution  equation  is  well  known.   If  K(t) 
is  the  input  waveform  to  a  linear  device  whose  impulse  response  is  h(t),  the  output  waveform  is  described 
using  the  convolution  equation, 


^(t)   =   fx(T)  h(t-T)dr  =J    x(t-T)h(T)  dx  (1-1) 


or  in  simplified  notation 

7(t)  =  x(t)  *  h(t)  <;i-2) 

where  "*"  denotes  convolution.   The  Laplace  transform  of  y(t)  is  given  by 

Y(s)   =  X(s)  H(s)  (^-3) 

which  shows  that  the  convolution  process  transforms  to  a  product  operation  in  the  s-domain.  The  decon- 
volution process  is  defined  as  the  solution  of  the  convolution  integral  equation  (1-1)  when  either  x(t) 
or  h(t)  is  unknown,  i.e.. 


h(t)   =  y(t)  (1/*)  x(t)  (1-^) 

H(s)   =   Y(s)  ^  X(s)  (1-5) 


t  Centre  de  Microelectronique  de  Grenoble,  Centre  National  d' Etudes  de  Telecommunications,  Meylan,  France. 


x(t)   =  y(t)  (1/*)  h(t)  (1-6) 

X(s)   =   Y(s)  T  H(s)  (1-7) 

Note  that  (1/*)  which  denotes  deconvolution  transforms  to  a  division  operation.   Both  cases,  i.e.,  the 
solution  for  h(t)  or  x(t) ,  are  basically  the  same  with  respect  to  the  deconvolution  problem. 

The  above  convolution  integral  relationship  is  defined  for  continuous  waveforms.   In  the  problems 
considered  here  the  waveform  data  are  not  continuous  but  are  discrete;  the  data  are  obtained  from  a 
sampling  process  followed  by  an  analog  to  digital  conversion. 

For  the  purposes  of  this  present  work,  only  discrete  variables  or  sequences  are  considered.   When 
the  question  of  error  is  considered,  it  is  relative  to  the  differences  between  one  sequence  and  another. 
The  error  due  to  representing  a  continuous  function  by  a  sequence  is  not  considered  here.   However,  it 
should  be  recognized  that  such  an  error  (aliasing)  can  be  very  significant  when  spectral  information  is 
desired  [1,2] . 

A  sequence  of  N-values  for  x,  x(0),  x(l) ,  x(2) ,  ...  x(N-l),  is  formally  denoted  as  {x(k)}  where  the 
k-th  number  in  the  sequence  is  x(k) , 

X  =   {x(k)}  (1-8) 

For  convenience,  here  the  notation  x(k)  will  denote  the  sequence  of  values  for  x. 

For  discrete  variables,  i.e.,  sequences,  the  convolution  and  deconvolution  equations  (1-1)  through 
(1-7)  are  respectively  replaced  by 


krl  k-1 

y(k)   =    2.  x(i)h(k-i)   =    Yl     x(k-i)h(i) 
i=o 


(1-9) 


i=o 


y(k)  =  x(k)*h(k)  (^_-L0) 

Y(s)  =  X(s)  H(s)  ^_-^^^ 

h(k)  =  y(k)  (1/*)  x(k)  (;L_3^2) 

H(s)  =  Y(s)  V  x(s)  (1-13) 

x(k)  =  y(k)(l/*)h(k)  ^_^^) 

X(s)  =  Y(s)  ^  H(s)  (-L_.L5) 

In  the  remainder  of  this  report,  for  the  most  part,  the  analyses  will  be  carried  out  in  terms  of  discrete 
variables.   Where  continuous  variables  are  used,  the  distinction  will  be  clearly  made.   However,  in  the 
figures,  for  convenience  in  drawing  the  discrete  waveforms  and  spectra,  they  will  be  drawn  as  continuous 
functions,  but  labeled  with  the  appropriate  discrete  variables. 


Numbers  in  brackets  indicate  literature  references  at  the  end  of  this  paper. 


In  the  experimental  situation  the  imperfections  of  the  acquisition  device  and  the  measurement 
system  noise  which  is  added  to  the  signal  make  the  available  data  not  representative  of  the  original 
continuous  waveform  but  of  something  close  to  it.   The  convolution  relationship  still  holds;  but,  in 
general,  the  deconvolution  may  no  longer  be  stable.   That  is  to  say,  the  solution  which  comes  from  the 
direct  brute-force  application  of  the  mathematical  equations  may  not  represent  the  desired  signal. 
Furthermore,  an  exact  solution  does  not  exist  and  must  be  replaced  by  an  estimate.   The  problem  is  to 
find  an  approximate  solution  which  is  stable  under  small  changes  in  the  initial  data.   To  do  that  one 
uses  a  regularizing  operator  or  filter.   This  filter  depends  on  the  characteristics  of  the  processed  wave- 
forms and  on  the  method  used  to  perform  the  deconvolution.   Particular  emphasis  is  given  to  the  fact  that 
the  experimenter  controls  the  operation  and  is  the  judge  of  the  quality  of  the  deconvolution  result. 

In  this  report  after  a  general  presentation  of  the  deconvolution  problem,  two  possible  solution 
methods  are  examined.   In  each  method  the  reasons  for  the  instability  are  inspected  and  lead  to  the 
synthesis  of  the  corresponding  regularizing  function.   The  experimental  results  include  the  evaluation 
of  the  quality  of  the  waveform  obtained  as  a  solution  of  the  deconvolution  process. 

Specifically,  this  report  is  divided  into  six  chapters,  the  first  being  the  present  introduction. 
Chapter  2  is  devoted  to  a  general  presentation  of  the  deconvolution  problem  and  possible  approaches  for 
solution. 

Chapter  3  describes  the  experimental  equipment  used  to  obtain  waveform  data  and  the  various  experi- 
mental quantitative  criteria  which  may  be  used  to  measure  the  quality  of  a  solution  to  a  deconvolution 
problem. 

Chapter  4  discusses  frequency  domain  deconvolution  methods  and  presents  detailed  examples  of 
solutions  to  deconvolution  problems  using  frequency  domain  processing  which  include  the  synthesis  of 
regularization  filters. 

Chapter  5  discusses  time  domain  deconvolution  methods  and  presents  several  examples,  using  time 
domain  processing,  one  illustrating  a  time  domain  solution  to  a  deconvolution  problem  which  was  also 
solved  in  Chapter  4  using  frequency  domain  methods.   The  synthesis  of  regularizing  filters  for  time 
domain  deconvolution  is  also  discussed. 

Chapter  6  presents  a  summary  of  the  deconvolution  methods  considered  and  some  general  conclusions 
about  the  solution  to  the  deconvolution  problem. 


2.   GENERAL  PRESENTATION  OF  THE  DECONVOLUTION  PROBLEM 
2.1   An  Exact  Solution  Exists  for  Exact  Waveforms 

Consider  the  linear  system  shown  in  figure  2-1;  if  the  input  waveform  x  (k)  and  the  measurement 
system  impulse  response  h  (k)  are  assumed  to  be  exactly  known,  then  the  exact  output  waveform  y  (k)  is 
also  known  and  is  given  by  the  convolution  equation, 

y^(k)   =  x^(k)  *   h^(k).  (2.1-1) 

Furthermore,  if  h  (k)  is  the  unknown,  the  deconvolution 

h^(k)   =   y^(k)  (1/*)  x^(k)  (2.1-2) 

is  a  unique  solution  which  exists  only  because  x  (k)  and  y  (k)  represent  exact,  error  free  waveforms 
(sequences) . 

2.2   Experimental  Situation  Yields  Inexact  Waveforms 

Practically,  the  exact  waveforms  x  (k)  and  y„(k)  are  not  available  from  measured  experimental  data, 
because  the  exact  values  are  corrupted  by  system  noise  and  sampling  errors.   Consequently,  the  experi- 
mental data  x(k)  and  y(k)  replace  the  exact  data  x  (k)  and  y  (k)  respectively.   However,  each  represen- 
tation, x(k)  or  y(k),  belongs  to  a  family  of  waveforms  defined  as  "close"  to  the  exact  waveform;  that  is 
to  say,  x(k)  belongs  to  the  family  [X], 

X  £  [X],  (2.2-1) 

such  that  for  all  elements  x  of  this  family, 

P(x,  X  )  <  6  (2.2-2) 

o   —   X 

where  p  being  defined  for  any  kind  of  waveforms  w  and  w  by  the  relation. 


'"1   ^  21 

P(w.%)  =  j^N   ^  ^"i  -  "oi^ 


1/2 
'  (2.2-3) 


The  parameter  6  may  be  interpreted  as  a  measure  of  the  error  which  is  made  when  the  waveform  w  is 
represented  by  w.   Then  for 

X  £  [X],  p(x,x^)  <  6  (2.2-4) 

y  £  [Y],  p(y,y^)  <  6^  (2.2-5) 

These  two  relationships  define  families  of  waveforms  such  that  "small  changes"  are  allowed  between  the 
elements  of  each  family. 

The  direct  use  of  the  deconvolution  equation  (1-3)  when  applied  to  the  y(k)  and  x(k)  waveforms 
leads  to  a  result  h(k)  which,  in  general,  is  not  close  to  h  (k)  the  exact  solution,  i.e.,  if  [H]  is  the 
family  of  waveforms  close  to  h  (k)  ,  then 

h(k)  i    [H]  (2.2-6) 

and 

p(h,  h^)>  6^  (2.2-7) 

hence,  the  solution  is  not  stable  under  small  changes  in  the  initial  data. 


A  strict  mathematical  solution  is  just  not  possible  because,  in  general,  the  result  is  not  stable  [3] 
The  problem  is  then  to  find  an  approximate  solution  stable  under  small  changes  in  the  initial  data.   The 
requirement  of  stability  has  to  be  satisfied  since  it  is  associated  with  the  physical  determinancy  of  the 
phenomenon  described  by  the  convolution  equation.   A  rule  has  to  be  chosen  for  the  definition  of  the 
family  [H]  of  possible  solutions.   In  practice  this  choice  is  arbitrary;  consequently,  one  needs  to  use 
supplementary  information  concernin  the  solution. 

However,  it  is  known  that  the  convolution  itself  is  stable  under  small  changes  in  the  initial  data. 
If  h(k)  is  the  obtained  approximate  solution,  its  convolution  with  the  known  input  waveform  x(k)  must 
give  a  result  which  belongs  to  the  family  [Y]  of  waveforms  close  to  y  (k) .   If 

c(k)  =  h(k)  *  x(k)  (2.2-8) 

and 

p(c,y)  <   6  ;  c  £  [Y]  (2.2-9) 

then 

e(k)  =  y(k)  -  c(k)  (2.2-10) 

defines  an  error  criterion  whose  properties  will  be  used  to  judge  the  quality  of  the  representation  of 
h  (k)  by  h(k).   Moreover,  it  may  be  assumed  that  the  approximate  solution  has  to  be  smooth  because  the 
solution  represents  a  continuous  deterministic  signal  or  system  (in  sampled  data  discrete  form.)   This 
qualitative  information  will  be  used  in  addition  to  the  error  criterion  to  control  the  process  which 
leads  to  the  accepted  solution. 


2.3  The  Regularization  Operator 

Another  way  which  employs  linear  algebra,  i.e.,  matrix  methods,  may  be  used  to  describe  the  convolu- 
tion process.   The  convolution 

y(k)  =  h(k)  *  x(k) 
may  be  written  as 

y(k)  =  A  x(k)  (2.3-1) 

where  A  is  an  operator  which  maps  the  family  [X]  into  the  family  of  waveforms  [Y] .   It  has  already  been 
pointed  out  that  the  inverse  operation  corresponding  to  the  operator  A   ,  which  was  expected  to  map  [Y] 
in  [X] ,  is  not  stable.   Hence,  A   has  to  be  replaced  by  a  new  operator  R  which  really  maps  [Y]  in  [X] 
such  that 

x(k)  =  R  y(k).  (2.3-2) 

R  is  called  the  "regularized  operator".   Its  expression  depends  upon  the  properties  of  the  two 
families  of  waveforms  [X]  and  [Y] .   These  properties  include  the  precision  in  the  knowledge  of  a  wave- 
form when  compared  to  the  exact  waveform  (parameter  6)  and  the  shape  of  the  waveform  itself.   On  the 
other  hand,  the  method  used  to  perform  the  deconvolution  determines  the  practical  effect  of  the  noise 
and  influences  the  definition  of  the  regularized  operator.   From  the  user's  point  of  view,  the  proper- 
ties of  R  must  be  easy  to  modify  in  such  a  way  that  a  progressive  improvement  of  the  regularization 
finally  leads  to  the  accepted  estimate  solution. 

Functionally, 

R  =  R  {[X],  [Y],  a(6)}  (2.3-3) 

where  a  is  a  set  of  parameters  which  are  used  to  take  into  account  the  different  possible  noise  figures. 


2.4  Possible  Methods  of  Deconvolution 

There  are  many  different  methods  for  solving  the  deconvolution  equation  [4,5].   In  this  section 
three  methods  are  briefly  discussed  in  order  to  provide  a  broader  perspective  of  deconvolution  methods. 
In  the  main  part  of  this  report  only  two  of  the  three  methods  will  be  discussed.   In  practice,  all 
methods  of  deconvolution  (including  the  three)  have  a  common  feature  in  that  the  amount  of  information 
that  can  be  extracted  from  the  empirical  curves  (physical  data)  is  limited.   Consequently,  some  form  of 
approximation  is  a  necessary  part  of  the  solution. 

2.4.1  Model  Fitting 

Model  fitting  consists  in  finding  an  analytical  expression  to  represent  the  exact  waveforms  y  (k) 
and  X  (k)  before  the  deconvolution  is  performed  e.g.,  Prony's  Method  [6,7].   In  this  case,  from  the 
inexact  available  data  it  is  necessary  to  obtain  a  mathematical  expression  which  is  assumed  to  represent 
the  exact  solution.   This  operation  has  to  be  made  for  each  waveform  with  the  additional  condition  that 
the  two  obtained  waveforms  have  to  be  compatible  in  the  convolution  process,  i.e.,  that  y'  (assumed 
exact  representation  of  y(k)  must  be  the  image  of  x'  (assumed  exact  representation  of  x(k))  in  the 
convolution  operation.   This  method,  model  fitting,  is  not  examined  in  this  present  experimental  study. 

2.4.2   Time  Domain  Processing 

The  problem  is  solved  using  only  the  time  domain  representation.   The  regularizing  operator  is  also 
defined  in  the  time  domain.   This  class  of  methods  includes  the  inversion  of  the  convolution  matrix  and 
the  iteration  technique. 

2.4.3   Transform  Domain  Processing 

Transform  domain  processing  methods  include  the  method  of  moments  (statistical)  [4],  the  homomorphic 
deconvolution  [8]  and  the  Fourier  transform  deconvolution.   The  latter  method  is  considered  here. 
Generally,  in  transform  domain  processing  the  available  data  are  transformed  into  another  representation 
(moments,  frequency,  etc.)  and  the  regularizing  operator  is  defined  in  this  new  domain  before  the 
return  to  the  time  domain. 

2.5  Loss  of  Information  in  the  Convolution  Process 

The  details  of  the  observed  waveform  are  smoothed  by  the  convolution  with  the  apparatus  function. 
It  is  impossible  to  distinguish  between  the  effects  of  the  noise  and  very  small  irregularities  due  to 
the  signal  itself.   Moreover,  the  amount  of  information  available  in  the  output  waveform  must  be  com- 
parable to  the  amount  of  information  needed  for  the  description  of  the  input  waveform  (or  impulse 
response) .   This  last  point  is  very  important  and  may  be  explained  using  a  frequency  domain  presentation. 

First  of  all,  when  the  duration  of  the  unknown  impulse  response  waveform  is  smaller  than  that  of  the 
input  waveform,  figure  2.2,  the  corresponding  discrete  frequency  domain  spectrum  magnitudes  are  as  shown 
in  figure  2.3.   It  is  evident  from  figure  2.3  that  the  discrete  spectra  X(n)  and  Y(n)  do  not  contain 
enough  high  frequency  components  to  allow  a  complete  description  of  H(n).   It  is  then  impossible  to 
determine  H(n)  or  h(t)  unless  additional  information  is  used. 


T'Jhen  the  duration  of  the  unknown  impulse  response  is  larger  than  that  of  the  input  pulse,  x(t)  and 
y(t)  contain  enough  information  to  determine  h(t),  figures  2-4  and  2-5.   Consequently,  a  possible  waj'  to 
remove  the  fundamental  difficulty  illustrated  by  figures  2-2  and  2-3,  is  to  use  an  input  signal  x(t)  of 
lesser  duration.   When  the  input  pulse  duration  is  less  than  that  of  the  system  impulse  response,  the 
only  problem  which  remains  is  to  find  a  way  to  manage  the  effects  of  the  noise  which  always  is  present. 
Such  management  is  called  a  regularizing  process  or  regularization. 

In  conclusion  of  this  chapter  on  the  general  presentation  of  the  deconvolution  problem,  it  is 
important  to  reiterate  (in  slightly  different  words)  a  point  made  in  an  earlier  study  [5] : 

"A  strict  mathematical  solution  cannot  be  obtained 
because  of  the  ever  present  noise  in  the  original 
data,  but  the  correlation  between  the  mathematics 
and  the  physical  data  must  alwa^'s  be  kept  in  mind 
to  reach  an  acceptable  estimated  solution". 

The  effect  on  the  final  solution  of  the  uncertainty  of  the  initial  data  due  to  noise  will  depend 
upon  the  method  used  for  the  deconvolution.   Two  methods  will  now  be  studied  in  detail:   the  Fourier 
transform  processing  method  and  the  time  domain  processing  method.   In  the  frequency  domain,  it  will  be 
seen  that  the  disturbing  effect  of  the  noise  appears  in  the  high  frequency  range.   In  the  second  method, 
i.e.,  in  the  time  domain  method,  it  will  be  seen  that  the  noise  increases  the  value  of  the  difference 
between  two  consecutive  samples,  and  this  random  deviation  from  the  exact  values  causes  instability  in 
the  solution.   Furthermore,  the  time  domain  method  is  very  sensitive  to  the  D.C.  value  or  offset  of  the 
waveform;  this  leads  to  difficult^'  in  using  the  method. 
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Figure  2.1   Linear  System  Possessing  an  Impulse  Response  h^(k). 
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Figure  2.2  The  unknown  h(k)  convolved  with  x(k)  where  the  duration  of 
h(k)  is  less  than  that  of  x(k). 


Figure  2.3  The  frequency  domain  spectrum  magnitudes  corresponding 
to  x(k),  h(k)  and  Y(k)  of  figure  2.2. 
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Figure  2.4   The  unknown  h(k)  convolved  with  x(k)  where  the  duration 
of  h(k)  is  greater  than  that  of  x(k). 
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Figure  2.5   The  frequency  domain  spectrum  magnitudes  corresponding  to 
x(k),  h(k)  and  y(k)  of  figure  2.4. 
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3.   EQUIPMENT  AND  METHODS  USED  IN  THE  EXPERIMENTS 
3.1  Acquisition  of  Waveforms 

A  minicomputer  controlled  sampling  oscilloscope  system  was  used  to  implement  the  deconvolution 
experiments.   The  specific  system  employed  was  the  NBS  Automatic  Pulse  Measurement  System  (APMS)  [9,10]. 
The  oscilloscope  bandwidth  is  very  large  (dc-18  GHz)  and  permits  the  observation  of  picosecond  domain 
signals.   A  block  diagram  of  the  APMS  is  shown  in  figure  3.1. 

In  operation,  a  repetitive  signal  synchronized  to  the  oscilloscope  time  base  is  applied  to  the 
sampling  head.   The  time  at  which  the  samples  are  taken  is  controlled  by  the  minicomputer  (using  the 
digital  to  analog  converter).   The  amplitude  of  each  sample  is  digitized  and  stored  in  the  memory.   The 
effect  of  the  noise  which  always  accompanies  the  signal  may  be  reduced  by  using  either  analog  or  digital 
averaging  techniques.   The  APMS  can  only  acquire  a  limited  finite  amount  of  data  (points);  consequently, 
a  time  window  of  finite  extent  is  used  to  observe  the  waveform.   The  acquisition  process  is  fast  and 
provides  a  large  number  of  data  points  (102A  for  example)  which  are  stored  in  an  array  whose  format  is 
well  adapted  for  further  digital  processing.   The  CRT  terminal  provides  a  powerful  tool  which  allows  the 
operator  to  implement,  control,  and  monitor  signal  acquisition,  conditioning  and  deconvolution  processing. 

3.2   Interactive  Processing 

Because  no  exact  solution  exists,  it  is  difficult  to  make  the  deconvolution  process  entirely 
automatic.   For  that,  one  would  have  to  use  a  general  mathematical  criterion  whose  definitions  could 
hardly  take  into  account  each  particular  configuration  of  available  data.   The  resulting  solution  would 
not  be  so  precise  as  the  one  which  is  obtained  if  the  judgment  of  the  operator  is  taken  as  a  part  of  the 
process.   With  the  operator  involved  in  the  process  the  graphical  capability  of  the  CRT  display  is 
extensively  used  to  maintain  the  contact  between  the  operator  and  the  digital  processing.   Each  time  a 
decision  is  made  its  effect  appears  in  the  display  permitting  a  measurement  of  its  effectiveness  or 
suitability. 

3.3  Qualitative  Measurement  of  the  Quality  of  One  Solution 

The  operator  uses  all  the  available  qualitative  information.   He  may  have  a  preconceived  idea  of 
the  shape  of  the  expected  waveform.   In  this  case  he  may  recognize  that  the  obtained  result  is  close  or 
far  off  the  solution  he  expects.   Often,  the  only  usable  criterion  is  the  smoothness  of  the  deconvolved 
waveform  which  results  from  the  deconvolution  process. 

3.4  Quantitative  Criteria 

Because  different  methods  may  be  used  for  the  achievement  of  the  deconvolution,  the  quality  of  the 
result  has  to  be  measured  using  a  quantitative  error  function.  This  criterion  permits  a  choice  between 
several  possible  solutions.   The  error  function  [e]  is  defined  by  the  relationship,  figure  3-2. 

[e]  =  [y]  -  {[d]  *  [x]}  (3.4-1) 

[y]  and  [x]  are  the  two  known  waveforms  (output  and  input),  and  [d]  is  the  approximate  solution  of  the 
deconvolution 

[d]  =  [y]  (1/*)  [x].  (3.4-2) 
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If  the  deconvolution  process  is  perfect  (exact  solution)  all  samples  of  the  error  function  [e]  will  be 
equal  to  zero. 

The  use  of  this  error  function  is  justified  by  the  following  remarks.   Remember  that  the  convolution 
operator  is  stable  under  small  changes  in  the  initial  data.   If  the  obtained  result  [d]  belongs  to  the 
family  [H]  of  waveforms  close  to  the  exact  solution  [h  ],  the  result  [y']  of  its  convolution  by  the  input 
waveform  must  belong  to  the  family  [Y]  of  waveforms  close  to  the  exact  output.   But  the  exact  output  [y  ] 
is  not  known.   Then  [Y]  is  used  as  the  reference  for  the  definition  of  the  error  function  (3.4-1). 

The  error  function  must  be  carefully  interpreted.   Four  parameters  will  be  used  for  its  characteri- 
zation, the  mean  value  e,  the  standard  deviation  a,    the  maximum  value  e    and  the  minimum  value  e  .  . 

max  mxn 

The  mean  value  is  given  by 

N 

1 

k=l 


J   E   e(k)  (3.4-3) 


Ideally  e  should  be  zero,  hence  in  practice  a  value  as  close  to  zero  as  possible  is  sought.   The  standard 
deviation  is  given  by 


N  _  - 

E   [e(k)  -  e]' 
k=l 


1/2 

(3.4-4) 


and  must  be  small  when  compared  to  the  maximum  value  of  the  reference  waveform  [y] .   A  typical  acceptable 
value  for  a   is  when  it  is  the  same  order  of  magnitude  as  the  noise  on  the  output  waveform  [y] .   Experi- 
mentally, it  has  been  shown  that  the  maximum  and  minimum  values  of  the  error,  e    and  e  .   occur  for 

max      mm 

some  samples  associated  with  the  largest  values  of  [y] ,  figure  3.2;   also,  for  a  well  chosen  regularizing 

function  le     and  e  .   will  be  of  the  same  order  of  magnitude.   Finally,  it  should  be  kept  in  mind 
'  max '      '  mm '  _ 

that  in  the  course  of  the  deconvolution  process,  the  goal  is  to  make  the  error  parameters  e,  a,  e    and 

f  '  1=  f  max 

e  .   very  small  compared  to  the  maximum  value  of  the  reference  waveform  [yl. 
mxn 
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Figure  3.1   The  NBS  Automatic  Pulse  Measurement  System 
(PTP  and  PTR  are  paper-tape-punch  and 
paper-tape-reader,  respectively) . 
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Figure   3.2     The  error   e(k)    is   the  difference  between  the  reference 
function  y(k)    and   the   estimated   function   d(k)    *  x(k). 
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4.   FREQUENCY  (TRANSFORM)  DOMAIN  DECONVOLUTION 

A.l   Introduction 

Frequency  domain  deconvolution  methods  operate  indirectly  on  the  time  domain  sequence  x(k) ;  they 
operate  on  the  transform  sequence  X(n)  which  is  the  Discrete  Fourier  Transformation  (DFT)  [11,12]  of  the 
time  domain  sequence.   The  DFT  and  its  inverse  IDFT  are  defined  in  the  following  way.   The  sequence  X(k), 
i.e.  , 

x(kT)   =  x(0),  x(T),  x(2T),  .  .  .  x[(N-l)T]  (A. 1-1) 

is  sequence  of  values  spaced  apart  by  an  interval  T.   The  DFT  of  x(k)  is  a  sequence  of  complex  values 
given  by 

N-1 
X(n  Q)   =    Yl       ^(kT)  e~^^  '^''",   n  =  0,  1,  2  ...  N-1  (A. 1-2) 

k=0 

where  each  value  is  spaced  apart  by  an  interval  Q,   where  fiT  =  2'n'/N.   The  inverse  transformation  IDFT  is 
given  by 

N-1 
x(kT)   =  ^   Z)  X(n  Q)    e^  ^  '^^'^ ,    k  =  0,  1,  2  ...N-1  (4.1-3) 

n=0 

As  indicated  above,  the  notation  x(k)  and  X(n)  implies  x(kT)  and  X(n  Q) ,    respectively. 
Errors  introduced  into  the  deconvolution  process  stem  from  three  sources: 

1.  In  the  time  domain  -  erroneous  or  inconsistent  data  and  computation 
errors  in  computing  the  error. 

2.  Transform  Domain  -  filtering  and  computation  errors. 

3.  Transformation  Operation  -  computation  errors. 

As  summarized  in  Chapter  2,  noise  or  errors  in  the  time  domain  data  cause  the  data  to  be  inconsis- 
tent in  the  sense  that  the  input  x(k)  and  output  y(k)  data  (i.e.  the  two  known  quantities)  are  not  re- 
lated to  each  other  through  a  causal  stable  impulse  response.   For  example,  if  for  a  given  causal  stable 
system  x  (k) ,  y  (k) ,  and  h  (k)  are  the  input,  output  and  system  impulse  response  sequences,  respectively, 
then  the  sequences  are  uniquely  related  by  the  convolution  equation  (1-8) . 

y  (k)  =  x  (k)  *  h  (k)  (4.1-4) 

o       o       o 

Consequently,  x  (k)  and  y  (k)  together  provide  the  data  to  exactly  characterize  the  system  h  (k) , 
figure  4.1.   Now  assume  that  in  the  measurement  of  x  (k)  and  y  (k)  the  independent  noise  sources,  n   and 
n  ,  corrupt  x  (k)  and  y  (k) ,  respectively,  to  yield  x(k)  and  y(k),  figure  4.2.   Clearly,  x(k)  and  y(k) 
are  different  from  x  (k)  and  y  (k) ,  respectively.   Furthermore,  in  general  x(k)  and  y(k)  are  not  related 
to  each  other  through  (4-4)  and  cannot  provide  h  (k)  through  the  deconvolution  process  as  specified  by 

h  (k)   =  y  (k)(l/*)  X  (k)  (4.1-5) 

o         o  o 

which  specifically  employs  x  (k)  and  y  (k)  and  not  x(k)  and  y(k) ,  respectively.   However,  the  impulse 

o         o 

response 

h(k)  =  y(k)(l/==--)  x(k)  (4.1-6) 

can  be  calculated.   This  result  may  be  stable  or  unstable,  and  if  stable,  it  may  resemble  h  (k)  or  a 
very  noisy  h  (k) . 
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The  impulse  response  of  a  physical  system  is  causal  (zero  for  time  less  than  zero,  the  time  at  which 
the  impulse  is  applied)  and  stable  (returns  to  zero  after  the  passage  of  time),  d  (t)  in  figure  4.3. 

By  operating  on  the  data  sequences  with  the  DFX  the  discrete  Fourier  transforms  of  x(k)  and  y(k)  can 
be  obtained;  these  transforms  are  the  complex  sequences, 

DFT[x(k)]  =  X(n)  (4.1-7) 

DFT[y(k)]  =  Y(n)  (4.1-8) 

Also,  it  is  true  that  the  exact  processes  (4.1-4)  and  (4.1-5)  have  the  DF  transforms 

Y  (n)  =  X  (n)  H  (n)  (4.1-9) 

o       o     o 

H  (n)  =  Y  (n)/X  (n)  (4.1-10) 

o       o     o 

but  the  ratio  Y(n)/X(n)  still  does  not  represent  H  (n) ,  i.e. 

o 

H  (n)  i   Y(n)/X(n)  (4.1-11) 


Furthermore, 


which  leads  to  the  conclusion 


H(n)  =  Y(n)/X(n)  (4.1-12) 


H(n)  i   H  (n)  (4.1-13) 

o 

Consequently,  the  DF  transformation  in  itself  does  not  remove  the  errors  that  change  X  (n)  and  Y  (n) 

into  X(n)  and  Y(n),  respectively. 

The  amount  of  error  in  erroneous  or  noise  corrupted  data  depends  upon  the  signal  to  noise  voltage 

ratio,  SNR,  e.g.,  [x(k)]    /a  ,  where  O  is  the  standard  deviation  (3.4-4)  of  the  zero-mean  measurement 
max  X         X 

system  noise  present  while  x(k)  was  being  measured.   For  a  very  large  SNR,  the  noise  has  a  very  small 
effect  upon  the  resultant  (sampled)  data  sequence.   With  decreasing  SNR's,  the  noise  becomes  increasingly 
dominant.   Typically,  for  the  most  part,  the  noise  spectrum  is  uniformly  distributed  (constant  spectrum 
amplitude)  and  is  weighted  or  filtered  by  the  network  properties  of  the  sampled-data  measurement  system. 
The  net  result  is  that  the  higher  frequency  components  of  a  low  pass  signal  well  within  the  measurement 
system  passband  are  increasingly  corrupted  by  noise  with  increasing  frequency.   Consequently,  the 
corrupted  DFT's  X(n)  and  Y(n),  respectively,  with  increasing  corruption  at  the  higher  frequencies. 
Furthermore,   their  ratio  gives  the  corrupted  transfer  fuction  H(n) ,  (4.1-12),  whose  IDFT  yields  the 
corrupted  impulse  response 


Tymi 

I  X(n)  I 


h(k)  =  IDFT  [H(n)]  =  IDFT^^  (4.1-14) 

I-     J 

which  may  or  may  not  be  stable,  and  which  may  not  even  resemble  the  exact  or  true  impulse  response 

sequence  h  (k) :   stable  and  unstable  responses  are  shown  in  figure  4.3. 
o 

-J  To  reduce  the  influence  of  noise  at  the  higher  frequencies  on  H(n)  and  ultimately  on  h(k),  H(n)  is 
filtered  with  a  low  pass  filter,  R(n),  which  is  judiciously  chosen  so  as  to  provide  a  stable  and  smooth 
impulse  response,  d(k).   R(n)  operates  on  H(n)  to  give  D(n), 

D(n)  =  R(n)  :  H(n)  (4.1-15) 

where  :  denotes  that  R(n)  operates  on  H(n) .   In  the  time  domain, 

d(k)  =  IDFT  D(n)  (4.1-16) 

The  error  e(k)  due  to  the  estimate  d(k)  is  then  obtained  from  (3.4-1)  and  characterized  as  described  in 
Section  3.4. 
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Too  much  filtering  will  yield  an  impulse  response  which  varies  too  slowly,  while  too  little 
filtering  may  yield  an  unstable  result,  figure  4.3.   The  choice  of  the  low  pass  filter  parameters  is  in 
fact  arbitrary  because  the  filtered  result  giving  the  impulse  response,  d(k),  is  subject  to  the  inter- 
pretation of  the  viewer.   Accordingly,  the  filtering  process  introduces  an  intrinsic  error  even  though 
filtering  may  be  necessary  to  achieve  a  stable  result,  i.e.,  a  causal  stable  impulse  response. 

In  addition  to  the  noise  corruption  of  the  time  domain  sequence  and  the  transform  domain  filtering 
distortion,  the  DFT  ,  IDFT,  and  other  computations  can  introduce  numerical  errors  which  also  affect  the 
choice  of  the  estimate  d(k).   The  DFT  and  IDFT  are  usually  implemented  using  one  of  the  fast  Fourier 
transform  (FFT)  algorithms  [13].   DFT,  FFT  and  other  operations  can  be  analyzed  to  estimate  the  resultant 
signal  to  noise  ratios  where  the  noise  is  due  to  numerical  truncation  and  rounding  required  by  finite 
computing  register  length  [14].   However,  the  errors  due  to  computations  are  relatively  small  compared 
to  those  due  to  noise  and  filtering. 

In  summary,  the  overall  deconvolution  processes  for  time  domain  sequences  using  transform  domain 
filtering  possesses  three  sources  of  error  whose  effects  can  only  be  reduced  but  not  eliminated.   Hence, 
it  is  clear  that  the  deconvolution  process  is  inexact,  i.e.,  the  true  answer,  h  (k)  cannot  be  obtained 
as  has  been  pointed  out  earlier  (Section  2.5),  and  reiterated  here. 

However,  because  "hope  springs  eternal  in  the  human  spirit"  techniques  have  been  and  continue  to  be 

developed  by  which  x(k)  and  y(k)  are  operated  upon  in  order  to  "determine"  h  (k) .   The  job  cannot  be 

done  in  an  exact  sense,  that  is,  h  (k)  cannot  be  determined,  but  it  can  be  estimated;  the  result  is  an 

o 

estimate  for  h  (k) .   In  many  practical  situations  the  results  are  good  enough  to  be  very  useful.   To  be 
good  enough  and  useful,  the  results  reasonably  approximate  some  preconceived  view  of  the  system  behavior 
derived  from  a  physically-based  mathematical  model.   The  argument  can  be  made  to  appear  more  plausible 
by  saying  that  the  effects  of  the  computation  noise  and  the  noise  components  n  (k)  and  n  (k) ,  figure  4.2 
are  relatively  small;  but  a  little  noise  can  ultimately  be  just  as  devastating  to  the  deconvolution 
process  as  a  large  amount  of  noise.   In  fact,  in  practice,  techniques  for  noise  reduction  or  enhancement 
of  the  signal  to  noise  ratio  are  first  used;  then,  that  noise  which  remains  will  still  be  potent  enough 
to  prevent  exact  deconvolution.   The  direct  effect  of  noise  can  be  seen  in  the  analysis  of  time  domain 
deconvolution  processes;  analyses  of  the  effect  are  presented  in  Chapter  5. 

4.2   Methods  Considered 

Two  frequency  domain  methods  are  considered  in  the  remainder  of  this  chapter;  they  are: 

1.  A  two  parameter  filtering  method 

2.  A  one  parameter  filtering  method. 

The  common  objective  of  the  two  methods  is  to  filter  out  the  noise  which  inhibits  the  deconvolution 
process  from  producing  a  smooth  stable  result.   The  filtering  process  is  called  a  regularizing  operation 
and  the  filter,  a  regularizing  filter,  is  denoted  by  R(n) . 

In  the  first  method  there  are  two  parameters  which  fix  the  characteristics  of  the  filter.   Its 
simple  mathematical  form  enables  the  viewer  to  easily  envision  and  choose  a  trial  set  of  values  for  the 
two  parameters  based  upon  spectrum  amplitude  data  in  dB,  20  log|H(n)|,  versus  frequency  n  Q.      After  the 
two  parameters  are  selected,  the  viewer  observes  the  time  domain  response  (IDFT),  d (k) .  of  the  filtered 
H(n) ,  D(n).   The  choice  of  the  filter  structure  in  this  method  was  not  based  upon  any  strong  a  priori 
conditions,  but  rather  upon  the  weak  condition  to  attenuate  the  noisy  portion  of  the  spectrum  for  H(n) 
and  thus  reduce  the  error  (3.4-1). 

In  the  second  method  there  is  only  one  filter  parameter;  its  value  governs  the  gradual  on-set  of 
filtering  to  provide  a  smooth  time  domain  response.   Again,  the  viewer  selects  a  trial  value  for  the 
filter  parameter  and  observes  the  time  domain  response,  d (k)  of  the  filtered  H(n) ,  D(n) .   However,  in 
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this  method  the  filter  structure  is  based  upon  the  a  priori  conditions  to  minimize  a  performance  measure 
consisting  of  a  weighted  sum  of  the  error  power  and  the  smoothness  power  (power  being  the  energy  per 
period  or  time  window.) 

4.3  The  Two  Parameter  Method 

4.3.1  The  Filter  Structure 

The  filter  specifically  operates  on  the  real  and  imaginary  parts  of  H(n) , 

H(n)  =  Hj^(n)  +  j  H^Cn)  (4.3-1) 

D(n)  =  Dj^(n)  +  J  D^(n)  (4.3-2) 


to  yield 


As  it  was  indicated  for  (4.1-14),  the  notation  R(n) :H(n)  denotes  that  R(n)  operates  on  H(n)  which  in 
general  implies  that  the  operation  need  not  be  the  product  R(n)H(n) .   In  the  two  parameter  method,  the 
operation  is  not  a  product  and  operates  on  H  (n)  and  H  (n)  to  yield 

K  1 


D,(n) 


H„(n),  0  <  n  <  n 
K       —       o 


|H(n  )|   exp(-Mn),   n  ^  n  £  N-1  (4.3-3a) 


H^(n)    ;   0  <  n  <  n 
I  —      o 

Dj(n)   =   I  (4. 3- 3b) 

o      ,   n   <  n  <  N-1 
o  —   — 


The  magnitude  of  D(n)  is  given  by 


[H^(n)  +  \i\{n)]^'^  ,      o  <  n  <  n^ 

|D(n) j  =   {  (4.3-4) 

lH(n  )|  exp(-Mn)  ,   n   <  n  <  N-1 
'o'    "^  o—   — 


and 


20  log[H^(n)  +  Hj(n)]^''^   ,   0  <  n  <  n^ 

D(n)|,^   =   {  (4.3-5) 

'  20  log|H(n  )1  -20(loge)Mn  ,   0  £  n  <  N-1 
I 
The  constant  M  is  the  cut-off  slope  of  filtered  H(n) ,  i.e.,  [D(n)i^g,  figure  4.4  and  is  given  by 

-100  +  iH(n  )  I 

-2_dB  (4.3-6) 


(A-l)n^ 

The  rationale  for  the  structure  of  the  filter  was  to  provide  a  simple  means  for  attenuating 
magnitude  lH(n)|  response  for  n  greater  than  some  value  n^  (i.e.,  in  the  noisy  region),  while 
simultaneously  totally  eliminating  the  noisy  phase  function,  e(n), 

.  Tan"^[H  fn)/H^(n)],   0  <  n  <  n 
)(n)   =  )        '     ^         -      o  (^.3_7) 

0  ,   n   <  n  <  N-1 
o  —   — 
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4.4  Application  of  the  Two  Parameter  Method 
4.4.1   Procedural  Summary 


In  practice,  the  viewer  inspects  the  data  H(n)    and  selects  the  acceptable  range,  n  ,  of  the  data 

do  o 

(relatively  noiseless  data);  this  specifies  n  and  |H(n  )|.   Then,  the  viewer  selects  the  cut-off  co- 
efficient, A,  which  then  completely  specifies  the  cut-off  slope,  M,  and  D(n).   Finally,  the  viewer  calls 
for  the  IDFT  of  D(n)  which  provides  the  estimated  or  deconvolved  impulse  response,  d(k).   If  the  resul- 
tant d(k)  and  its  associated  error  criteria,  (3.4-1)  through  (3.4-4),  are  not  satisfactory  in  the  judg- 
ment of  the  viewer,  the  process  is  repeated  using  different  values  of  n  and  A. 

The  procedure  starting  with  the  time  domain  sequences  x(k)  and  y(k)  and  ending  with  the  deconvolved 
time  domain  impulse  sequence  d(k)  progresses  as  follows: 

1.  Offset  correction  for  x(k)  and  y(k).   Eliminate  the  offset,  i.e.,  insure 
that  x(o)  =  x(N-l)  =  0  and  y(o)  =  y(N-l)  =  0.   If  x(k)  and  y(k)  are  step- 
like waveforms  so  that  x(N-l)  and  y(N-l)  are  not  zero,  they  must  be  trans- 
formed to  impulsive  waveforms.   See  subsection  4.6.2. 

2.  DF  transform  x(k)  and  y(k)  to  obtain  X(n)  and  Y(n);  display  |x(n)|    and 

dB 

l^(-)ldB- 

3.  Calculate  the  complex  ratio  H(n)  =  Y(n)/X(n);  display  |H(n)|   . 

an 

4.  Select  the  regularization  filter  coefficients  n  and  A;  then  filter 

o 

H(n)  to  yield  D(n),  (4.3-3)  display  |D(n)L^. 

dB 

5.  IDF  transform  D(n)  to  obtain  the  time  domain  deconvolved  result, 
d(k).   (3.4-2),  (4.1-16);  display  d(k). 

6.  Calculate  the  error  criteria  (3.4-1),  (3.4-3),  and  (3.4-4)  display  e(k) 

and  the  values  for  e    ,  e  .   ,  e,  and  a. 
max .   mm . 

7.  Judge  the  quality  of  the  deconvolved  impulse  response  estimate,  d(k), 
in  terms  of  the  error  criteria;  if  unsatisfactory,  return  to  step  4 
and  change  the  filter  parameters. 

8.  If  the  quality  is  satisfactory,  print  the  d(k)  graphic  display, 
related  error  criteria,  and  the  numerical  values. 

Examples  of  the  two  parameter  deconvolution  method  are  given  below  for  the  determination  of  the  impulse 
response  of  (1)  300  meter  coaxial  transmission  line  and  (2)  a  very  broadband  antenna. 

4.4.2   The  Estimation  of  the  Insertion  Impulse  Response  of  300  Meters 
of  RG  58  C/U  Coaxial  Cable 

In  this  example  the  insertion  input  and  output  waveforms  of  a  300  meter  length  of  RG  58  CU  coaxial 
transmission  line,  x(k)  and  y(k),  respectively,  were  recorded,  figure  4.5;   the  signal  to  noise  ratios 
for  these  waveforms  were  relatively  large,  SNR-x  75  dB  and  SNR-y  46  dB.   Deconvolution  to  obtain  the  in- 
sertion impulse  response  was  then  performed  using  varying  degrees  of  regularization  with  the  two  parameter 
filter.   The  DFT  of  the  insertion  impulse  response  is  the  scattering  insertion  parameter  S^_(n).  [17]. 
The  discrete  transformation  and  inverse  transformation  operations  were  implemented  using  the  FFT  and 
IFFT  algorithms. 

Next,  various  combinations  of  SNR-x  and  SNR-y  were  created  by  adding  noise  (a  pseudorandom 
sequence)  into  the  original  high  SNR  data  recorded  above.   For  each  combination  of  SNR's,  two  parameter 
deconvolution  was  performed,  again  using  varying  degrees  of  regularization. 
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Specifically,  the  following  deconvolution  experiments  were  performed: 

1.  SNR-x:   75  dB,  SNR-y :  46  dB 

No  regularization:  (N-1)  Q   =  250  MHz 

Moderate  regularization  =  n  fi  =  100  MHz,  A  =  5 

Strong  regularization:   n  Q  =  74  MHz,  A  =  6 

o    o  o 

2.  SNR-x:  75  dB,  SNR-y:  36  dB 

No  regularization:  (N-1)  fi  =   250  MHz 
Strong  regularization:   n  fi  =  74  MHz,  A  =  6 

3.  SNR-x:   45  dB,  SNR-y:  46  dB 

No  regularization:  (N-1)  Q  =  250  MHz 

Moderate  regularization:   n  Q  =   100  MHz,  A  =  5 

4.  SNR-x:  40  dB,  SNR-y:  40  dB 

No  regularization:  (N-1)  Q  =   250  MHz 

Strong  regularization:   n  Q  =  74  MHz,  A  =  6 

The  results  are  summarized  in  figures  4.6  through  4.8  and  more  detailed  information  is  presented  in 

figures  4.9  through  4.23. 

Figure  4.6  brings  together  the  results  for  the  no-regularization  cases  of  each  of  the  experiments. 

In  each  of  these  cases,  H(n)  was  directly  inversely  transformed  by  the  IDFT  to  obtain  the  time  domain 

sequence  h(k).   Also,  shown  in  figure  4.6  are  the  |H(n)|    corresponding  to  the  various  h(k)  responses. 

do 

Generally,  it  is  seen  that  the  noisier  a  |H(n)|    is,  the  noisier  is  the  corresponding  h(k).   Further- 

dB 

more,  a  change  in  the  SNR-y  of  the  output  waveform  is  more  critical  than  that  of  the  input  waveform  with 

regards  to  the  noise  characteristics  of  |H(n)|    and  h(k).   For  example,  a  10  dB  reduction  in  SNR-y 

do 

yields  very  noisy  results,  figure  4.6b,  while  a  30  dB  reduction  in  SNR-x  yields  a  much  smaller  noise 

increase  in  |H(n)|    and  h(k),  figure  4.5c.   This  follows  from  the  fact  that  in  the  latter  case  the 
dB 

transmission  properties  of  H(n)  attenuates  the  effect  of  the  noise  components  introduced  at  the  input. 
Figure  4.7  brings  together  the  results  for  the  regularized  cases  of  each  of  the  experiments.   In 
each  of  these  cases,  H(n)  is  operated  on  by  the  two  parameter  (n  ,A)  regularization  operator  to  yield 
D(n),  and  the  result  is  inversely  transformed  by  the  IDFT  to  obtain  the  time  domain  sequence,  d(k). 
The  horizontal  scale  is  expanded  to  50  ns/Div.  (as  compared  to  200  ns/Div.  in  figure  4.6)  in  order  to 
show  the  fine  structure  of  d(k).   Also,  h(k)  for  the  largest  SNR  case  (figure  4.5a)  is  shown  in  ex- 
panded time  scale  as  figure  4.5a  for  comparison  purposes. 

It  is  evident  that  the  regularization  or  filtering  process  produces  much  smoother  estimates,  d(k), 
of  the  actual  coaxial  cable  impulse  response,  h  (k) ,  than  do  the  various  non-regularized  h(k).   The 
smoothest  non-regularized  response  is  shown  in  figure  4.7a  while  the  other  non-regularized  responses 
of  lesser  smoothness  are  shown  in  figure  4.6.   The  effect  of  regularization  is  very  dramatic  in 
reducing  the  noise  from  the  impulse  response. 

Before  commenting  any  further  on  figure  4.7,  it  is  necessary  to  consider  the  input  and  output 
waveforms  of  figures  4.9  and  4.10,  respectively.   First,  from  figure  4.9,  it  is  clear  that  the  input 
pulse  waveform  abruptly  rises,  at  say  t^ ,  and  is  zero  prior  to  t^;   also,  the  pulse  duration  was  of  the 
order  of  10  ns .   Second,  from  figure  4.10,  the  output  pulse  waveform  possesses  some  small  variations 
prior  to  the  large  abrupt  rise  of  the  pulse  waveform;  also,  there  is  a  small  variation  from  the  smooth 
decay  of  the  pulse  waveform  at  about  800  ns.   The  main  rise  and  subsequent  return  to  zero  of  the  out- 
put waveform  occurs  in  about  1000  ns ;  consequently  the  input  waveform  is  approxmiately  an  impulse, 
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which  in  turn,  says  that  the  output  waveform  will  approximate  the  impulse  response  of  the  coaxial  cable, 

h  (k) .   Equipped  with  these  physical  details  of  the  input  and  output  waveforms,  further  comments  can  be 
o 

made  about  the  deconvolution  results  presented  in  figures  A. 6  and  A. 7. 

Because  the  output  waveform  in  figure  4.10  was  produced  by  a  very  short  duration  input  pulse, 
variations  similar  to  those  in  the  output  waveform  should  be  detectable  in  the  deconvolved  waveforms, 
d(k)  and  h(k).   Inspection  of  figure  4.7  shows  different  oscillations  before  the  main  rise  of  each 
waveform.   The  time  scale  in  figure  4.7  is  50  ns/Div.  which  is  one-quarter  of  that  in  figure  4.10. 
Because  figure  4.10  doesn't  possess  oscillations  in  the  50  ns  time  Interval  prior  to  the  start  of  the 
main  waveform  rise,  those  appearing  in  figure  4.7  are  due  to  noise  or  the  filtering  process. 

For  example,  refer  to  figure  4.6a  and  note  the  noise  level  on  the  baseline  of  the  waveform;  then 
look  at  figure  4.7a  which  is  an  expanded  time  scale  version  of  figure  4.6a.   Consequently,  it  is  clear 
that  the  oscillations  before  the  main  rise  are  due  to  noise. 

Next,  look  at  figure  4.7b  and  note  that  the  regularization  has  smoothed  out  the  variations  due  to 
noise,  but  that  the  smoothing  very  close  to  the  initial  rise  is  not  as  great  as  that  elsewhere. 
Stronger  regularization  increases  the  magnitude  of  the  oscillations  as  evidenced  in  figure  4.7c.   This 
is  due  to  the  filtering  by  the  regularization  process  which  introduces  the  Gibbs  phenomena  due  to  the 
truncation  of  H(n)  when  forming  D(n)  [15].   Such  effects  can  be  made  negligible  by  avoiding  abrupt  or 
rectangular  filtering  of  the  H(n)  spectrum.   Smooth  filtering  is  used  in  the  one-parameter  deconvolution 
method  discussed  in  Section  4.6. 

In  Section  4.1  it  was  pointed  out  that  the  FFT  and  IFFT  computations  temselves  introduce  noise  or 
errors  due  to  finite  register  length.   Such  errors  also  cause  the  imaginary  part  of  the  IFFT  to  be 
nonzero.   For  example,  say,  the  real  time  domain  sequence  f(k)  is  given  and  is  then  transformed  by  the 
FFT  to  F(n);  next,  f(k)  is  recovered  from  F(n)  via  the  IFFT.   However,  upon  inspection  of  the  data  for 
IFFT  {F(n)}  it  is  found  that  the  result  is  complex  with  a  very  small  imaginary  part  while  the  real  part 
is  very  close  to  the  values  of  the  original  sequence  f (k) .   In  practice,  such  non-zero  imaginary  parts 
are  always  present,  and  are  relatively  very  small,  figure  4.8.   In  Section  4.5,  the  nature  and  origin  of 
the  errors  affecting  d(k)  will  be  discussed. 

4.4.3  The  Estimation  of  the  Insertion  Impulse 
Response  of  a  Broadband  Antenna 

As  a  second  example  for  the  application  of  the  two  parameter  method,  the  method  was  applied  to  ex- 
perimental data  to  obtain  the  insertion  impulse  response  of  a  broadband  antenna.   The  data  was  obtained 
using  the  NBS  Time  Domain  Antenna  Range  in  conjunction  with  the  NBS  APMS.   Schematically,  the  experimen- 
tal arrangement  is  shown  in  figure  4.24. 

Figures  4.25  and  4.26  show  the  input,  x(k),  and  the  output,  y(k),  and  the  magnitudes  of  their 
respective  DFT's,  jx(n)|  and  |Y(n)|.   The  magnitude, 

|H(n)|    =   |Y(n)|  ^  |x(n)| 

is  shown  in  figure  4.27.   Deconvolution  results  are  shown  in  figures  4.28a  and  4.28b  for  two  different 
degrees  of  filtering;  note  that  the  smoother  result  and  most  error  occurs  for  the  greater  degree  of 
filtering. 
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4.5  On  the  Nature  of  the  Error,  e(k) 

With  the  results  presented  in  Section  4.4  in  mind,  it  is  now  appropriate  to  provide  a  model  for  the 
error  e(k)  encountered  in  frequency  domain  deconvolution.   The  model  is  not  restricted  to  any  specific 
frequency  domain  regularization  filter;  consequently,  it  is  valid  for  both  the  two-parameter  method 
(Sections  4.3,  4.4)  and  the  one  parameter  method  (Section  4.6).   The  discussion  will  be  divided  into 
three  parts: 

1.  Error  in  the  absence  of  both  regularization  (filtering)  and  computation 
errors  (noise) . 

2.  Error  in  the  presence  of  regularization  and  absence  of  computation  errors. 

3.  Error  in  the  presence  of  both  regularization  and  computation  errors, 
being  contained  in  subsections  4.5.1  through  4.5.3,  respectively. 

4.5.1   Error  in  the  Absence  of  the  Filtering  Operation  and  Computation  Noise 

The  error  sequence,  e(k),  as  given  by  (3.4-1)  includes  the  effects  of  filtering  and  computation 
errors  in  the  d(k)  term, 

e(k)  =  y(k)  -  x(k)  *  d(k)  (4.5.1-1) 

In  the  absence  of  filtering  and  computation  noise 

d(k)  =  h(k)  (4.5.1-2) 

then 

e(k)  =  0  (4.5.1-3) 

because 

y(k)  =  x(k)  *  h(k)  (4.5.1-4) 

In  words,  the  error  sequence,  e(k),  will  be  zero  only  when  d(k)  equals  h(k);  any  other  choice  for  d(k) 
will  give  a  non-zero  result. 

4.5.2   Error  in  the  Presence  of  the  Filtering  Operation 
and  in  the  Absence  of  Computation  Errors 

Now  consider  the  nature  of  h(k)  for  the  signals  x(k)  and  y(k)  which  are  both  noisy  signals,  i.e., 

x(k)  =  X  (k)  +  n  (k)  (4.5.2-1) 

o        X 

y(k)  =  y^(k)  +  n  (k)  (4.5.2-2) 

where  n  (k)  and  n  (k)  are  the  noise  components  associated  with  the  true  signals  x  (k)  and  y  (k) , 

X         y  o         o 

respectively.   The  system  impulse  response  is 

h(k)  =  [y^(k)  +  n  (k)](l/*)[x^(k)  +  n^(k)]  (4.5.2-3) 


whose  DFT  is 


Y  (n)  +  N  (n) 

O  X 
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The  transfer  function  H(n),  (4.5.2-4),  obtained  from  the  deconvolution  y (k) (l/*)x(k) )  can  be 
expressed  as 

1  +  N  (n)/Y  (n) 

""^^^    -     "o^-)  1  ^  N\n)/X°(n)-  (A.5.2-5) 

X      o 

where 

Y^(n) 

^o^"^  =  T~M  (4.5.2-6) 

Define,  ° 


1  +  N  (n)/Y  (n) 
^^"^  =    1  +  S^(n)/X^(n)  (4.5.2-7) 

H  (n)  is  the  DFT  of  the  stable  true  impulse  response  defined  by 

h^(k)  =  y^(k)(l/*)x^(k)  (4.5.2-8) 

Consequently,  H(n)  can  be  viewed  as  the  product  of  the  true  transfer  function  H  (n)  and  some  noisy 
transfer  function,  M(n) ,  which  may  or  may  not  be  stable. 

H(n)  =  H  (n)  M(n)  (4.5.2-9) 

o 

If  H(n)  is  filtered  by  R(n)  to  deemphasize  the  effect  of  the  noise,  M(n) ,  then  D(n)  results 

D(n)  =  R(n)  [H  (n)  M(n)]  (4.5.2-10) 

o 

Assuming  that  R(n)  operates  only  on  M(n) ,  i.e.,  write 

D(n)  =  H  (n)  [R(n):M(n)]  (4.5.2-11) 

o 

in  the  limit 

lim    D(n)  =  H  (n)  (4.5.2-12) 

R:M^  1         ° 

However,  it  is  physically  impossible  for  the  filter  R(n)  to  operate  solely  on  M(n).  In  practice  R(n) 

also  operates  on  H  (n)  with  the  design  of  R(n)  being  such  that  it  affects  M(n)  more  than  H  (n)  and 
makes  d(k)  at  least  stable  and  possibly  smooth. 

Consequently,  filtering  attempts  to  make  h(k)  stable  by  remvoing  m(k)  from  the  impulse  response, 

h(k)  =  h  (k)  *   m(k)  (4.5.2-13) 

o 

in  such  a  way  that 

d(k)  =  h  (k)  *  [r(k)  *   m(k)]  (4.5.2-14) 

o 

that  is,  d(k)  becomes  approximately 

d(k)  =  h  (k)  (4.5.2-15) 

o 
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where  the  last  equation  is  a  judgment  by  the  individual  performing  the  deconvolution.   The  resultant  d(k) 
does  not  cause  the  error,  e(k),  to  go  to  zero;  only  h(k)  can  do  so. 

Therefore,  it  can  be  concluded  that  a  paradoxial  situation  exists:   h(k)  makes  the  error,  e(k), 
vanish;  but  h(k)  is  not  the  desired  result,  h  (k) ,  which  the  individual  performing  the  deconvolution 
wishes  to  achieve.   The  paradox  occurs  because  y(k)  is  used  as  the  reference  waveform  or  sequence 
in  the  error  expression,  (3.4-1).   As  pointed  out  in  Section  3.4,  y  (k)  is  not  known;  consequently, 
the  known  quantity  y(k)  is  used. 

Since  h(k)  and  H(n)  are  given  by  (4.5.2-3)  and  (4.5.2-4),  respectively,  the  error,  e(k),  can  be 
expressed  as 

e(k)  =  y(k)  -  x(k)*r(k)  *  [y (k) (l/*)x(k) ]  (4.5.2-16) 

=  y(k)  *  [6(k)  -  r(k)]  (4.5.2-17) 

where  6(k)  is  a  unit  impulse.   In  the  discrete  frequency  domain,  the  error  is 


E(n)  =  [Y^(n)  +  Ny(n)]  -  R(n) 


rY  (n)  +  N  (n)-] 

X  (n)  +  N'(n)    f^o(->  ^  \^-^^ 
o  X    J 

(4.5.2-18) 


which  corresponds  to  (4.5.2-16)  while 


E(n)  =  [1  -  R(n)]  [Y  (n)  +  N  (n) ]  (4.5.2-19) 

corresponds  to  (4.5.2-17).   The  square  of  the  error  magnitude  is 

|E(n) 1^  =   il  -  R(n) 1^  [y  (n)  +  N  (n) 1^  (4.5.2-20) 

and  the  error  energy  over  the  N  points  comprising  the  time  window  or  one  period  of  e(k)  is,  by  the 
Parseval  relation  [22,  23], 


N-1  N-1 

E   e^(k)  =     ^        Y.       |E(n)r  (4.5.2-21) 


where 


N 
k=0  n=0 


N-1  N-1 

y      |E(n)!    =     y        !l  -  R(n)  r  |y  (n)  +  N  (n) r  (4.5.2-22) 

n=0  n=0  ^ 


Consequently,  it  is  seen  that  the  error  energy  per  period  depends  upon  the  filter  characteristics  R(n) 

Y  (n) ,  and  N  (n) .   The  \ 
o  y 

of  (3.4-4),  is  given  by 


Y  (n) ,  and  N  (n) .   The  variance  or  the  square  of  the  standard  deviation  of  the  error,  i.e.,  the  square 
o  y 


,      T    N-1  _ 

O       =     ^       Y,       [e(k)  -  e  ]^  (4.5.2-23) 

k=0 
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?      1^-1?       -7 
a''  =  i   y;   e  (k)  -  e   "■  (4.5.2-24) 

2     1     ^"1   ,     ,2-2 
a^  =  ^    ^    |E(n)r   -  e  (4.5.2-25) 

N     n-0 

where  e  is  the  mean  value  of  the  error  (3.4-3). 

The  effects  of  filtering  on  the  error  can  be  graphically  illustrated  by  considering  the  magnitude 
of  the  error,  E(n).   From  (4.5.2-20)  the  magnitude  is 

E(n) I   =   ll  -  R(n) I   |y  (n)  +  N  (n) I  (4.5.2-26) 

Figure  4.29  shows  the  error  for  two  choices  for  R(n)  where  R(n)  is  real;  the  absolute  error,  |E(n)I, 
approaches  zero  as  R(n)  approaches  unity. 


4.5.3   Error  in  the  Presence  of  the  Filtering  Operation  and  Computation  Errors 

In  the  presence  of  filtering  and  computation  errors  the  estimated  impulse  response,  d(k),  in  (3.4-1) 
is  replaced  by  d.(k).   The  error  is  then  given  by 

e(k)  =  y(k)  -  x(k)  *  d2(k)  (4.5.3-1) 

which  can  be  explicitly  written  as 

e(k)  =  y(k)  *    [6(k)  -  r(k)  -  m^(k)]  -  x(k)  *  [r(k)  -  n^^(k)  +  n^  (k)  +  n^^{]!.)]    +   n^(k) 

(4.5.3-2) 

and  which  will  be  derived  below.   e(k)  reduces  to  (4.5.2-13)  when  the  computation  errors  vanish,  that  is, 
when 

m  (k)  =   6(k),  unit  impulse  (4.5.3-3) 

nj^^(k)   =  0  (4.5.3-4) 

n^^Ck)   =  0  (4.5.3-5) 

n  (k)    =0  (4.5.3-6) 

e 

The  derivation  of  e(k),  (4.5.3-2),  and  the  definitions  of  (4.5.3-3)  through  (4.5.3-6)  are  based 
upon  the  assumption  that  the  round-off  error  in  calculating  a  specific  quantity,  can  be  represented 
as  an  additive  term  to  the  desired  quantity  [16].   For  example,  if  A(n)  is  being  calculated  from  a(k) 
by  the  FFT,  then  the  result  is 

A^(n)   =   A(n)  +  N^(n)  (4.5.3-7) 

where   A(n)  is  the  desired  result,  and  N'(n)  is  the  computation  error  or  noise.   Similarly,  the  IFFT 
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of  A, (n)  would  also  introduce  an  additive  error,  n  . (k) . 
1  al 


a^(k)   =   a(k)  +  n  (k)  +  n  ^ (k) 

1  a       al 


(4.5.3-8) 


Thus,  it  is  assumed  that  the  cummulative  error  depends  upon  the  number  of  computations,  and  that,  a(k) 
is  not  exactly  recoverable  from  the  numerical  calculations. 

To  derive  e„(k),  (4.5.3-2),  the  starting  point  is  x(k)  and  y(k)  and  the  procedure  follows  below. 


FFT  [x(k)]   =   X  (n)   =   X(n)  +  N' (n) 


(4.5.3-9) 


FFT  [y(k)]   =  Y^(n)   =  Y(n)  +  N^(n) 


(4.5.3-10) 


Here,  N' (n)  and  N' (n)  are  computation  errors  (noise):  note  that  both  X(n)  and  Y(n)  also  contain  measure- 

X         y 

ment  noise  corresponding  to  the  noise  contained  in  x(k)  and  y(k),  (4.5.1-5)  and  (4.5.1-6),  respectively. 
Next  compute  the  quotient  of  (4.5.3-10)  divided  by  (4.5.3-9);  this  result  is  defined  as  H  (n) , 


H^(n) 


Y^(n) 


+  KAn) 


X^(n)      hi' 


(4.5.3-11) 


where  N,  ^ (n)  is  the  additive  error  due  to  the  calculation  of  H. (n) . 
hi  1 

Filtering  H  (n)  yields  D  (n) . 


D^(n)   =   R(n) 


Y^  (n) 


+  N^i (n) 


X  (n)      hi 


+  N  (n) 


(4.5.3-12) 


where  N  (n)  is  the  additive  error  due  to  the  calculation  of  D.  (n) .   In  terms  of  the  theoretical  trans- 
forms Y(n)  and  X(n),  D  (n)  is  given  by 


D^(n)   =   R( 


,rY(n)  j  1  ^  N;(n)/Y(n)  ) 
"^[^X(n)  i  1  +  N^(n)/X(n)  ^'  ^  ^■ 


(n) 


+  N^(n) 


(4.5.3-13) 


^("^[IS)-  >^i(->^\i(")]  ^\(-) 


(4.5.3-14) 


where  M  (n)  is  a  multiplicative  error  term  defined  by  the  {    ;  term  in  (4.5.3-13). 
The  IFFT  of  D  (n)  is  given  by 


d2(k)  =  r(k)  ■"|[y(k)(l/*)  x(k)]*  m^(k)  +  ti^-^'^''^  }'^   "r  ^''^  "^  ''d2  ^^^ 


(4.5.3-15) 


where  n   (k)  is  a  complex  additive  error  due  to  the  IFFT  calculation  of  d  (k) .   The  error  e(k)  is  given 
by 


and 


e(k)  =  y(k)  -  x(k)  *  d2(k)  +  n^(k) 


d_(k)  =  d,(k)  +  n  (k)(l/*)x(k) 
5  I  e 


(4.5.3-16) 
(4.5.3-17) 


then 


e(k)  =  y(k)  -  x(k)  *  d2(k) 


(4.5.3-18) 
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e(k)  =  y(k)  *  [6(k)  -  r(k)  *  m^ (k) ]  -  x(k)  *  [r(k)  *  n,  . (k)  +  n^(k)  +  n,_(k)]  +  n  (k) . 

X  n±       r       dz        e 

(4.5.3-19) 

where  n  (k)  is  the  additive  error  due  to  the  calculation  of  e  (k)  . 
e 

When  the  filtering  operation  is  deleted 

r(k)  =   6(k)  (4.5.3-20) 

n  (k)  =  0  (4.5.3-21) 

d_(k)  will  change  to  d,(k)  and  the  error  becomes 

e(k)  =  y(k)  -  x(k)  *  d^(k)  (4.5.3-22) 

Explicitly,  from  (4.5.3-19)  the  error  becomes 

e(k)  =  y(k)  *  [6(k)  -  m^(k)]  -  x(k)  *  l%^(^)   +  n^2(^)l  +  n^(^)-  (4.5.3-23) 

This  is  the  error  obtained  when  no  filtering  is  employed  and  computation  noise  is  present.   When 

computation  noise  is  absent,  (4.5.3-3)  through  (4.5.3-6)  prevail  and  e(k)  reduces  to  zero. 

As  pointed  out  in  connection  with  (4.5.3-15),  d„(k)  is  complex  due  to  the  complex  additive  error 

term,  n,.(k),  which  results  from  the  IFFT  calculation.   Consequently,  the  estimated  impulse  response, 
dz 

(4.5.3-15) 

d2(k)  =  r(k)  *  [y(k)(l/*)x(k)  *   m^(k)  +  nj^^(k)]  +  n^(k)  +  n^^*^'^) 

or 

d2(k)  =  r(k)  *  [h(k)  *  m^(k)  +  ti^3^(k)]  +  n^(k)  +  n^^^^^)  (4.5.3-24) 

is  complex.   Also,  the  estimated  impulse  response,  without  filtering,  d,(k),  is  complex.   From 
(4.5.3-24)  with  (4.5.3-20)  and  (4.5.3-21)  d  (k)  is  given  by 

d^(k)  =   [h(k)  *  m^(k)  +  n^^W]    +   n^(k)  +  n^^  (k)  (4.5.3-25) 
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4.6   The  One  Parameter  Method 
4.6.1  The  Filter  Structure 

This  filter  is  an  optimum  filter  in  the  sense  that  it  minimizes  a  performance  measure  with  respect 
to  the  weighted  sum  of  the  error  power  and  a  smoothness  power  [18,19].   It  turns  out  that  there  is  only 
one  adjustable  parameter  in  this  method,  and  the  parameter  adjusts  the  weighting  between  the  error  and 
smoothness.   The  derivation  of  the  filter  is  given  below. 

Given  the  two  signals  x(k)  and  y(k) , 

x(k)   =  X  (k)  +  n  (k)                                 (4.6.1-1) 

y(k)   =  y^(k)  +  n  (k)                                 (4.6.1-2) 
computation  and  filtering  lead  to  the  estimated  impulse  response,  d„(k),  (4.5.3-16), 

d2(k)   =   r(k)  *  [y(k)(l/*)  x(k)  *   m^(k)  +n^^(k)]  +  n^ (k)  +n^2(^)               (4.6.1-3) 
and  the  error  function. 

e(k)   =  y(k)  -  x(k)  *  d2(k)                            (4.6.1-4) 

In  addition  to  the  error  constraint,  e (k) ,  define  an  auxiliary  constraint,  s(k), 

s(k)   =  c(k)  *  d2(k)                                  (4.6.1-5) 

where  c(k)  is  the  discrete  function  that  defines  the  constraint.   In  the  present  discrete  case  c(k)  is  the 

second-difference  operator,  which  corresponds  to  the  second  derivative  operator  for  continuous  functions, 

2    2 
d  /dt  ,  thus  s(k)  depends  upon  the  rate-of-change  of  the  first  derivative  of  d. (k) ,  i.e.,  its  smoothness. 


N-1 
E  =     XI  [e(k)]  ,  (4.6.1-6) 

k=0 


_      N-1 

S  =    E  [s(k)]''  (4.6.1-7) 

k=0 

Now,  to  obtain  an  optimum  solution,  adjust  d  (k)  until  the  performance  measure,  P, 

P  =  E  +  Y  s"  (4.6.1-8) 

is  minimized.   Y  is  a  weighting  function  which  selects  the  degree  of  smoothing.   Consequently,  the 
minimization  condition  is  approximately  expressed  by 


3d2(k) 


E  +  Y  S  I 


(4.6.1-9) 
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assuming  that  N  is  a  very  large  number  so  that  d-(k)  approximates  a  continuous  function  rather  than  a 
discrete  one. 

In  the  DF  domain  d  (k) ,  e (k) ,  and  s(k)  become 


°2(")   =   ^(">[lS 


j   M^(n)  +   N^^(n) 


E(n)   =  Y(n)  -  X(n)D2(n) 


S(n)   =   C(n)  D2(n) 


+  N  (n)  +  N^^(n) 
r       d^ 


(4.6.1-10) 
(4.6.1-11) 
(4.6.1-12) 


Parseval's  relation  [22,23] 


E  a  (k) 
k=l 


1   N-1 
n=0 


(4.6.1-13) 


jives  for  E  and  S, 


N-1 


n=0 


E(n)  E  (n) 


(4.5.1-14) 


N-1 


V  S(n)  S"(n) 


n=0 


(4.5.1-15) 


where  E  (n)  denotes  the  complex  conjugate  of  E(n),  etc.   The  sum  E  +  Y  S  is  given  by 

N-1 


^  +  "^  ^   "   N    ^     ^^^^^    E  (n)  +  Y  S(n)  S  (n)  ] 


n-0 


(4.5.1-16) 


where 


E(n)  E  (n)   =   [Y(n)  -  X(n)  D2(n)]  [Y(n)  -  X(n)  D2 (n) ] 


(4.6.1-17) 


S(n)  S  (n)   =   C(n)  D2(n)  C  (n)  D2(n) 


|C(n)  D2(n) 


(4.6.1-18) 


Consider  E(n)  E  (n) ; 


E(n)  E  (n)   = 


[Y(n)  -  X(n)  D2(n)]  [Y(n)  -  X(n)  D2(n)] 


[Y(n)  -  X(n)  X)^M]     [Y"(n)  -  {X(n)  D2(n)}"] 


=   Y(n)  Y"(n)  -  Y(n)  {X(n)D^(n)}"  -  Y"(n)  {X(n)D  (n)} 


+   {X(n)  D2(n)}  {X(n)  D2(n)} 
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E(n)E*(n)  =  |Y(n)|^  -  Y(n)  {X(n)  D2(n)}*   -  Y*n)  {X(n)  D^Cn)}   +   |x(n)  D(n)|^ 

(4.6.1-19) 
Hence 


E  +  Y  S   =  I     2  I  |Y(n)  1^  +  |x(n)  D^Cn)  |^  +  Yic(n)  D^Cn)  |^  -  Y^di)  {X(n)  D^di)} 
n-0  ' 

-  Y(n)  {X(n)  D^Cn)}*!  (4.6.1-20) 

Consider  X(n)  D  (n) .   From,  (4.6.1-10), 

X(n)  D^dn)   =   R(n)  [Y(n)  M^(n)  +  X(n)  Nj^^(n)]  +  X(n)  [N^(n)  +N^2(n)]         (4.6.1-21) 

First,  Im. (n) I  is  of  the  form  1  +  J(n),  where  J(n)  is  complex  and  |j(n)|<<l.   Second,  |n,  ,(n)l,  In  (n) I , 
1  hi        r 

and  |N,„(n)|  are  the  magnitudes  of  computation  errors;  consequently,  they  too  will  be  much  less  than 
unity.   When  X(n)  D„(n)  is  written  in  the  form 

X(n)  D2(n)   =   R(n)  Y(n)  A(n)  (4.6.1-22) 

A(n)  will  be  complex  and  of  the  form 

A(n)   =   1  +  K(n)  (4.6.1-23) 

where 

|K(n)|  «   1  .  (4.6.1-24) 

Consequently,  A(n)  will  be  very  close  to  the  real  number  unity. 
Some  complex  relations: 

XD  =  RYA  (4.6.1-25) 

=    IrIIyIIaI       ■*     +4     +  d)  (4.6.1-26) 

.     r  y  a 

(XD)*     =    IrI  IyI  IaI    /-((t>     +4)       +  <t>   )  (4.6.1-27) 

^      ^  I     I  I     I  I     1   /       ^r  y __a 

Y     =      |y|    /4)  (4.6.1-28) 


Y 


YI/-4  (4.6.1-29) 


■^1 


y""'(XD)    =    IrMyI     lAJ/cf)     +4)  (4.6.1-30) 

(XD)*  =  1r1  |y1^|a|  J-{<p^  + 


Y    (XD)*  =    IrMyI^IaI    /-((})     +  c})    )  (4.6.1-31) 


Y*(XD)  +  Y(XD)""'   =  |r1|y|^|a|  {2  Cos  (<[)   -  <f  ) )  (4.6.1-32) 

r    a 

=  2R|y|^!a|  Cos  (j)  (4.6.1-33) 


where 


Putting  (4.6.1-33)  into  (4.6.1-20)  yields 


(4.6.1-34) 
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N-1 
E+YS      =     I        X!     ||Y(n)|^   +    |R(n)    Y(n)    A(n)I^   +  Y|C(n)    D^CnXl^ 


n=0 


-2R(n)    Y(n)  |      A(n)      Cos    (})( 


»>} 


(4.5.1-35) 


i     I;'    |Y(n)|2-.    |R(n)Y(n)A(n)i2+Y|^(^i^y^(^^^-A(nl 
n=o 


-  2R(n)   Y(n)        A(n)    Cos   0( 


»') 


(4.6.1-36) 


i        E      I^WlM    l^R^n)     |A(n)|^+^    |C(n)|^R^(n)     |A(n)|^ 

"      --'^  '  |x(n)r 


n=0 


2   R(n)     |A(n)  ]    Cos   <t>in}\ 


(4.6.1-37) 


In  principle,  the  performance  measure  (4.6.1-8) 

P  =  E  +  Y  S" 

can  be  expressed  in  terms  of  discrete  time  domain  quantities  or  in  terms  of  discrete  frequency  domain 
quantities  by  virtue  of  the  Parseval  relation  (4.5.2-21).   However,  in  practice,  given  the  time  domain 
quantities  x(k)  and  y(k) ,  subsequent  numerical  computations  introduce  errors  into  the  frequency  domain 
quantities.   Consequently,  in  the  frequency  domain  the  result  is  ?', 

P'  =  E'  +  Y  ^' 

which  is  not  equal  to  P.   To  minimize  P,  d  (k)  is  adjusted.   Varying  d  (k)  corresponds  to  adjusting  the 
frequency  domain  filter  R(n) .   Consequently,  P  and  P'  are  minimized  through  the  relations, 

9  P        ^ 


d2(k) 


3  P' 


a  R(n) 


=  0 


respectively.   Because  the  computation  errors  are  very  small,  a  first  order  approximation  is  to  assume 
that  P  is  ec 
R(n),  i.e.. 


that  P  is  equal  to  P'.   With  that  assumption  being  made  a  change  in  d.(k)  corresponds  to  a  change  in 


and 


3d2(k) 


[^  +  ^0=  3fao["^^'] 


-,  ^,  ,  i  E  +  Y  S  I  =  .i    I,   .]        E  +  Y  sl  +  1^-7-       E  +  Y  S 
9  R(n)  {_     'J     3 1  R(n)  1  L        J     3  <t>_  L        J 


(4.6.1-38) 


(4.5.1-39) 


Consequently, 
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C^[e  +  Ys]   =  i   E    |Y(n)|2  (2iR(n)i!A(n)|2   +   2y  |  C(n)A(n)  |  ^  JRCn) 


n=0  |X(n) 

-  2|A(n)  |Cos[(})^(n)  +  (f)^(n)]  +  2  |  A(n)R(n)  |  Sin[(t)^(n)  +  (J)^  (n)  ]  i  (A. 6. 1-40) 

To  minimize  [E  +  Y  S]  with  respect  to  R(n)  ,  lR(n)|  and  (f)  (n)  must  be  such  that  (4.6.1-40)  is  zero. 

I     1 2 
Since  |Y(n)]   is  not  generally  zero,  the  {    }  -term  in  (4.6.1-40)  must  be  made  equal  to  zero  to  force 

(4.6.1-40)  to  zero.   Choosing  the  filter  phase  angle  to  be  the  negative  of  the  computation  error  phase 

angle  <t>^M  . 

())   (n)   =  -<p    (n)  (4.6.1-41) 

eliminates  the  phase  angle  dependence,  and  gives  for  |r  (n) | , 

|R^(n)|   =  ^^^ T  (4.6.1-42) 

lA(n) I |x(n) I   +  Y|A(n) | |c(n) \ 

Consequently,  the  optimum  filter  |r  (n)  |   /(})   (n)  is  given  by  (4.6.1-41)  and  (4.6.1-42)  in  the  presence 
of  computation  noise. 

If  smoothing  is  not  employed,  7  is  set  equal  to  zero;  and  the  filter  becomes 


|R  (n)  I     =   u/.  I  (4.6.1-43) 

'  o    'y=o      |A(n) I 


[*ro^")lY=0  ^  *a^"^  (4.6.1-44) 

which  says  that  the  filter  magnitude  and  phase  are  such  that  they  just  cancel  out  the  computation  errors 
so  as  to  force  the  error  function  e(k)  to  zero.   Furthermore,  when  computation  errors  are  not  present 
A(n)  is  unity,  i.e., 

|A(n)|   =  1  (4.6.1-45) 

cf)    =0  (4.6.1-46) 

which,  in  turn,  means  that  (4.6.1-43)  and  (4.6.1-44)  reduce  to 


[R(n)]   =  1  (4.6.1-47) 

Y=0 
A(n)=l 


[<})„(n)]  =  0  (4.6.1-48) 

^  Y=0 


A(n)=l 

which  is  the  case  of  no  filtering,  i.e.,  R(n)  =  1.   As  R(n)  approaches  unity,  the  filtering  l-R(n) 
approaches  zero,  figure  4.29. 
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In  the  presence  of  smoothing,  Y  ^  0,  and  in  the  absence  of  the  computation  errors,  A(n)  =  1,  the 
optimum  filter  is  real  and  is  given  by 

I     1 2 

R(n)   =  ^^^^ (A. 6. 1-49) 

|x(n)r  +  7|C(n)r 

Before  discussing  the  form  of  C(n),  some  summary  comments  are  in  order  to  clarify  the  meaning  of  the 
optimal  result  in  the  presence  of  computation  errors,  (4.6.1-41)  and  (4.6.1-42).   It  must  be  kept  in 
mind  that  in  any  given  problem  A(n)  is  not  known,  i.e.,  |A(n)|  and  (\>    (n)  are  unknown;  consequently,  the 
optimal  filter  cannot  be  calculated.   Equations  (4.6.1-41)  and  (4.6.1-42)  are  only  of  theoretical 
interest  in  that  they  show  that  computation  errors  require  that  the  optimum  filter,  R  (n) ,  be  different 
from  that  in  the  absence  of  computation  errors,  R(n) .   Specifically,  R  (n)  equals  A(n)  R(n) .   In  practice, 
since  A(n)  is  not  known,  it  is  convenient  to  let  R(n)  be  real  and  given  by  (4.6.1-49),  which  is  the  case 
in  the  absence  of  computation  errors.   Again,  because  A(n)  is  very  close  to  the  real  value  of  unity,  the 
use  of  (4.6.1-49)  is  the  best  that  can  be  done. 

c(k)  is  defined  as  the  second  difference  sequence  operator  and  corresponds  to  the  second  derivative 
operator  in  the  continuous  case.   The  first  backward-difference  sequence,  V,  is  defined  by  the  operation 
{x(kT)  -  x[k-l)T]}/T  where  the  interval,  T,  between  the  k-th  and  (k-l)-th  values  is  explicitly  shown. 
Implicitly,  the  operation  is  denoted  as  x(k)  -  x(k-l).   V  *  x(k)  is  given  by 


V  *  x(k)  =  [1,-1,0,...,  N-1]  *  x(k)  (4.6.1-50) 

2 
The  second  backward-difference  sequence,  V  ,  is  defined  by  the  implicit  operation  x(k)-2x(k-l)+x(k-2) . 

2 
V  "'  x(k)  is  given  by 

V^  *  x(k)   =   [1,-2,1,0,...,  N-1]  (4.6.1-51) 

2 
Note  that  there  are  N-3  and  N-4  zeros  in  the  V  and  V   sequences,  respectively.   The  DFT  of 

V^  =   [1,  -2,  1,  0,  0,  0,  N-1.]  (4.6.1-52) 


is  given  by  the  complex  sequence 


A^  =   ]^   V^(k)e  ^"^'^'^,  n  =  0,1,2,0,0,0,...N-1  (4.6.1-53) 


k=0 


1  -  2e"^^  '^^  +  e"J^  "^"^  ;   n  =  0, 1,  2,0,0, 0,  ..  .N-1  (4.6.1-54) 


which  also  contains  N-4  zeros,  and  where 


n  T  =  -i^  (4.6.1-55) 

N 


2 
Since   c(k)    is   defined   as   V    (k) ,    then   C(n)    is    given   by 


.    2   TTn  .    4   TTn 


C(n)      =      l-2e  +     e  ,    n=   0, 1 , 2 ,0,0,0, . . .N-1  (4.5.1-56) 
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and 


C(n)         =   C(n)C(n) 


.  2    TTn  .4    Tin 

,         "^    N        _^     "-'    N 
l-2e  +  e 


.  2   TTn  .4    TTn 

,    .    ^    N  ,         ^    N 

l-2e  +     e 


n  =   1,1,2,0,0,0,.. .N-1 


(4.5.1-57) 


|c(n)|^      =      6     -8   Cos  ^-^     +     2   Cos  ^Y^   ;      n  =  0,1, 2,0,0,0, ..  .N-1 


which  can  be  reduced  to  [21] 


.  4 


C(n)    =   Sin  (2TTn/N);   n  =  0,1, 2 , 0,0,0, .  .  .  ,N-1 


(4.6.1-58) 


(4.6.1-59) 


When  N  becomes  very  large,  the  argument,  2TTn/N,  approaches  a  continuous  function  of  f;   also,  the 
argument  becomes  very  small,  i.e.. 


large 

N 


L   N   J 


2  TTf,   0  <f<  1. 


(4.6.1-60) 


Since 


Sin  oj  =  oj-  —^  + 


3    5    7 

0)    .  OJ      w 


'      7  ' 


1\ 


+  •  •  •  ;  oj  <  CO 


(4.6.1-61) 


then  for 

^^^^^  [C(n)|^  =   |c(jw)|^  =  0)^  (4.6.1-62) 

which  corresponds  to  the  square  of  the  magnitude  of  the  continuous  function  second  derivative  operator. 


2        2 
D^  =  (j  oj)^ 


(4.6.1-63) 


dt 


D    =  w 


(4.6.1-64) 
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4.6.2  Conversion  of  a  Step-Like  Waveform  into  a  Duration  Limited  One 

Given  a  waveform  or  sequence,  f(k),  comprised  of  N  values,  successive  application  of  the  DFT  and 
IDFT  yield  an  F(j(jj)  and  f  (k)  which  are  periodic  of  period,  N.   Two  examples  are  shown  in  figure  4-30, 
where  the  given  sequence  of  N  points  are  at  the  left-hand  side  of  the  figure.   The  resultant  IDFT's  are 
periodic  in  which  one  cycle  is  identical  to  the  original  f(k). 

For  the  sequence  f  (k) ,  the  spectrum  F  (n)  corresponds  to  the  intended-time  sequence  f  (k) .   In 

the  case  of  f,  (k)  the  spectrum  G  (n)  does  not  correspond  to  the  intended-time  sequence.   Consider  the 
b  b 

following  physical  situation,  the  step-like  waveform  f,  (k)  is  N  points  of  a  sequence  that  begins  at 

zero  and  continues  to  n  =  >»,  remaining  at  the  constant  value  B  for  k  J^  k.  where  k  £  N-1.   Because  the 

DFT  will  use  only  N  points,  the  f,  (k)  is  abruptly  truncated  at  k  =  N-1.   The  resultant  spectrum  G,  (n)  is 

b  b 

not  the  intended-spectrum,  F  (n) .   The  IDFT  shows  g^ (k) ,  a  cycle  of  which  is  not  the  intended  f,  (k) , 

b  b  b 

i.e.,  the  constant  step-like  waveform.   The  abrupt  truncation  f.  (k)  introduces  the  spectral  components 

b 

of  an  abrupt  transition  (step  function). 

A  step-like  sequence,  f(k),  which  has  achieved  a  constant  value,  B,  for  k.  £  N-1,  can  be  converted  to 
a  duration-limited  one,  of  length  M  =  2N,  figure  4.31.   The  conversion  is  accomplished  by  extending  and 
delaying  f (k)  by  N  units  to  obtain  f(k-N),  figure  4.31a  and  b.   f(k)  and  f(k-N)  both  consist  of  M  =  2N 
points.   f(k-N)  is  subtracted  from  f (k)  to  yield  f„(k). 

f  (k)   =   f(k)-f(k-N),   k  =  0,1,2, •••M  -1  (4.6.2-1) 


The  DFT  of  f„(k)  is  given  by 


F2(m)   =   F(m)  -  F(m)e  J^^™^,  m  =  0,1,2, •••,  M-1       (4.6.2-2) 


2tt 
where   F(m)    is    the   DFT   of    f(k)    and    Q  =  —   .      F    (m)    is    the   product 


where 


Because 


F^Cm)    =    (1  -   e  i^'^^^)    F(m)  m  =   0,1,  2  ,  •  •  •  ,M-1  (4.6.2-3) 

F    (m)    =   F    (m)    F(m)  (4.6.2-4) 


F^(m)    =   1   -  e-^"^"^  =   1  -   e-^'"^  m  =   0,1,2,  •••,   M-1  (4.6.2-5) 


|f    (m)|      =      2    |Sin  mTT/2|  (4.6.2-6) 


L(m)      =      Tan-1       TT^^  (4.6.2-7) 

1  1-Cos   mTT 


It^/nI  0,m=0,2,4--*  /■/coQ\ 

I F    (m)         =  (4.6. 2-0) 


1 


2,    m  =   1,3,5- 


|)^(m)      «     0,   m  =   0,1,2, •••  (4.6.2-9) 
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the  DFT's  of  F(n)  and  F  (m)  are  simply  related.   Figure  4-32  shows  F(n)  and  F  (m)  where  m  has  been  used 
in  the  argument  of  F  .   Note  that  [ F  (m) j  is  twice  the  value  of  |F(n)|  for  m  =  2n-l.   Explicitly, 
F(n)  and  F  (m)  are  F(ni;)  and  F  [m(n/2) ] .   The  discrete  frequency  nfi  corresponding  to  n  =  N-1  is  equal  to 
the  discrete  frequency  mfi/2  with  m  =  M-2  because  M  =  2N.   Thus,  the  spectrum  of  F  (m)  has  (a)  the  same 
number  of  non-zero  values,  N,  as  that  for  F(n),  but  are  graphically  spaced  twice  as  far  apart,  and  (b) 
the'magnitude  of  JF.(m)[  is  twice  that  of  |F(nn)|  at  the  non-zero  values. 

In  deconvolution  studies  the  effect  of  the  time  domain  conversion  of  step-like  waveforms  for  x(k) 
and  y(k)  yields  an  estimated  impulse  response  due  to  a  unit  impulse  applied  at  k=0,  followed  by  a  nega- 
tive unit  impulse  applied  at  k=N-l,  figure  A. 33. 
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A. 6. 3  Application  of  the  One  Parameter  Method 
4.6.3.1  Procedural  Sunimary 


In  the  discussion  presented  here,  the  notation  will  be  that  employed  in  Section  4.5;  consequently, 
d(k;)  will  be  more  precisely  specified  by  d„(k);  similarly  D^ (n)  will  replace  D(n),  and  so  on  for  other 
quantities  as  the  occasion  arises.   The  reason  for  doing  so,  is  to  emphasize  that  the  quantities  contain 
errors.   In  the  sections  before  4.5,  the  errors  were  not  delineated. 

The  procedure  starting  with  the  time  domain  sequences  x(k)  and  y(k)  and  ending  with  the  deconvolved 
time  domain  impulse  sequence  d  (k)  progresses  as  follows: 

1.  Offset  correction  for  x(k)  and  yCk).   Eliminate  the  offset,  i.e.,  insure 
that  x(o)  =  x(N-l)  =  0  and  y(o)  =  y(N-l)  =  0.   The  offset  error  is 
systematically  introduced  by  the  data  acquisition  device.   If  x(k)  and  y(k) 
are  step-like  waveforms  so  that  x(N-l)  and  y(N-l)  are  not  zero,  they  must 
be  transformed  to  impulsive  waveforms.   See  subsection  4.6.2. 

2.  DF  transform  x(k)  and  y(k)  to  obtain  X, (n)  and  Y  (n) ;  display  [x  (n)l    and 

11  1     do 

|Y^(n)|^g,  (4.5.3-9),  (4.5.3-10). 

3.  Calculate  the  complex  ratio  H  fn)  =  Y^(n)/X^(n),  (4.5.3-11);  display  |h  fn)i  ,^. 

Ill  1    dB 

4.  Select  the  regularization  filter  coefficient  y;    then  filter  H  (n)  to  yield 
D^(n),  (4.5.3-14)  display  |D^(n)|^g. 

5.  IDF  transform  D  (n)  to  obtain  the  time  domain  deconvolved  result, 

d2(k).   (4.5.3-14);  display  d  (k) . 

5.   Calculate  the  error  criteria  (4.5.13-19);  display  e(k)  and  the  values  for  e     , 
_  max. 

e  .   ,  e,  and  0. 
mm. 

7.  Judge  the  quality  of  the  deconvolved  impulse  response  estimate,  d  (k) ,  in 
terms  of  the  error  criteria;  if  unsatisfactory,  return  to  step  4  and  change 
the  filter  parameter. 

8.  If  the  quality  is  satisfactory,  print  the  d  (k)  graphic  display,  related  error 
criteria,  and  the  numerical  values. 

Three  examples  of  the  one-parameter  deconvolution  method  will  now  be  presented;  they  are: 

1.  The  estimation  of  the  insertion  impulse  response  of  300  meters  of  RG  58  C/U  coaxial 
cable;  this  example  uses  the  same  insertion  waveform  data  as  was  used  in  the  two- 
parameter  method  example.  Section  4.4.1. 

2.  The  insertion  impulse  response  of  a  filter  for  a  step-like,  100  ps  transition 
duration,  reference  waveform  generator. 

The  first  example  allows  the  reader  to  compare  the  one  and  two  parameter  methods;  the  second 
example,  in  addition  to  illustrating  the  one  parameter  method,  illustrates  the  procedure  for  using 
step-like  waveform  data  rather  than  duration  limited  waveforms. 
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4.6.3.2  Estimation  of  the  Insertion  Impulse  Response 
of  300  Meters  of  RG  58  C/U  Coaxial  Cable 

The  insertion  waveforms  x(k)  and  y(k)  are  the  same  ones  used  in  the  two  parameter  method  experiments 
in  Section  4.4.2,  figures  4.9  and  4.10,  respectively.   Three  degrees  of  regularization  were  used,  Y  =  100, 
500,  and  10,000.   The  estimated  insertion  impulse  responses  and  their  corresponding  error  functions,  e(k), 
and  filtered  (regularized)  spectrum  magnitudes  are  shown  in  figures  4.34  through  4.36.   The  regularization 
filter  was  that  of  (4.6.1-49), 

„,  ,         |X(n) 1^ 

R(n)  =  '  %'  ' r 

lx(n)r  +  Y|c(n)r 

which  is  a  real  filter. 

When  Y  =  100,  the  regularization  is  moderate.   Compare  figure  4.34  with  figures  4.7b  and  4.12  of 
the  two  parameter  method.   The  results  are  comparable.   Increasing  Y  to  500  provides  strong  regulariza- 
tion.  Compare  figure  4.35  with  figures  4.7c  and  4.13  of  the  two  parameter  method;  note  that  the  Gibbs 
phenomena  has  appeared  in  the  results  of  the  two  parameter  method,  figure  4.7c,  while  the  one  parameter 
method  produces  no  Gibbs  oscillations.   Increasing  Y  to  10,000  provides  still  stronger  regularization 
with  no  evidence  of  Gibbs  oscillations  (compare  figures  4.36  to  4.7c  and  4.13).   Furthermore,  figure 
4.36  very  clearly  shows  the  small  variations  (reflections)  discussed  in  the  sixth  paragraph  of  section 
4.4.2.   Thus,  it  is  seen  that  the  one  parameter  method  uses  a  smooth  filtering  that  does  not  introduce 
the  Gibbs  effect  (as  was  pointed  out  earlier  in  section  4.4.2). 


38 


4.6.3.3    The  Insertion  Impulse  Response  of  a  Filter  for  a  Step-like, 
100  ps  Transition  Duration,  Reference  Waveform  Generator 

In  this  example,  the  insertion  impulse  response  of  a  liquid  filled  coaxial  transmission  line  is 
estimated.   The  coaxial  line  is  about  20  cm  long  and  is  filled  with  a  dilute  binary  solution  in  which 
heptance  is  the  solvent  and  heptanone  is  the  solute.   Because  heptanone  Is  a  polar  molecule  (heptane 
is  non-polar),  the  solution  comprises  a  dilute  polar  dielectric  and  as  such,  exhibits  a  dielectric 
relaxation  law.   Consequently,  the  frequency  dependence  of  the  dielectric  constant,  e(s)  causes  the 
transmission  line  to  behave  as  a  low  pass  filter.   The  filter  is  used  to  shape  a  20  ps  transition 
duration  (10%-90%)  step-like  waveform  into  a  100  ps  transition  duration  waveform  [25]. 

The  insertion  step-like  waveforms  x(k)  and  y(k)  are  converted  to  duration  limited  waveforms,  X-.  (k) 
and  y^ (k) ,  and  the  DFT's  X  (n)  and  Y  (n)  are  obtained,  respectively,  figure  4.37.   The  effects  of  filter- 
ing on  the  DFT's  are  shown  in  figure  4.38;  here,  the  real  filter,  R(n),  of  (4.6.1-49), 

Rfn^   -         |X(n)|^ 
R(n)   -  ■ y' r 

ix(n)|^  +  7|c(n)r 

2    4       6  —  — 

has  been  used  for  7  =  0,  10  ,  10  ,  and  10  .   R(n)  operates  on  H  (n)  to  yield  D  (n) 

D^(n)  =  R(n)  H^(n) 
and  the  IDFT  of  D  (n)  is  d„(k).   The  estimated  impulse  response  for  very  strong  regularization,  Y  =  10 


is  shown  in  figure  4.39.   Note  that  d  (k) 


d2(k)  =  y(k)(l/*)  x(k) 


is  the  response  due  to  a  pair  of  impulses  as  shown  in  figure  4.33,  and  d  (k)  is  the  desired  estimated 
impulse  response,  the  first  pulse  of  the  pair  in  d„(k). 

This  example  completes  the  discussion  of  frequency  domain  deconvolution.   In  the  next  chapter, 
some  time  domain  deconvolution  methods  will  be  presented. 
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K^(k) 


y„(k) 


Figure  4.1 


X  (k)  and  y  (k)  are  the  exact  sequences  which  are 
uniquely  related  through  the  system  impulse 
response  sequence  h  (k) . 


x^(k) 


Measurement 
System 


t 

x(k) 


h_,(k) 

1 

1 

n^(k) 


y„(k) 


Measurement 
System 


t 

y(k) 


n^(k) 


Figure  4.2 


x(k)  and  y(k)  are  corrupted  versions  of  x  (k)  and  y  (k) 

The  noise  sources  n  (k)  and  n  (k)  are  independent  of 

X         y 

each  other. 
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d^(k) 


d2(k) 


d3(k) 


d^(k) 


^*  k 


|-^-  k 


n-1 


Figure  4.3  The  effects  of  filtering,   d  (k) ,  the  stable  response  resulting 
from  the  "proper"  amount  of  filtering,  h  (k)  in  the  operator's 
judgment,   d  (k) ,  almost  stable,  needs  more  filtering.   d^Ck), 
very  unstable,  very  little  filtering,  if  any.   d  (k) ,  stable, 
too  much  filtering,  response  is  slowed  down. 


41 


0  dB 


H(no) 


100  dB 


o  o 

Discrete  Frequency,  nf.; 


Figure  4.4  The  magnitude,  |D(n)|j_,,  of  the  filtered  H(n)  versus  n  Q, 

dB 

in  dB,  20  log  |D(n) | ,  and  the  magnitude  of  the  unfiltered 
|H(n)1^3. 
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Figure  4.5   Insertion  response  measurement  method.   The  observables  are   the 

waveforms  x(k)  and  y(k)  obtained  before  and  after  the  system  under 
test  is  inserted  between  the  50  9,   impedance  generator  and 
measurement  system,  respectively.   H(n)  =  Y(n)/X(n)  =  S_  ,  the 
insertion  scattering  parameter  [17]. 
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Figure  4.6 


Deconvolution  without  any  regularization,  i.e.,  full  spectrum  of  H(n)  is 

used;  h(k)  is  shown  for  various  input  and  output  signal  to  noise  ratios. 

(a)  SNR-x  75  dB,  SNR-y  46  dB.   (b)  SNR-x  75  dB,  SNR-y  36  dB.   (c)  SNR-x 

45  dB,  SNR-y  46  dB.   (d)  SNR-x  40  dB,  SNR-y  40  dB.   |H(n)L^  is  shown  in 

db 

the  upper  right  hand  portion  of  each  figure  with  the  vertical  scale  of 
-10  dB/division  commencing  at  0.0  dB  and  the  horizontal  scale  of 
25  MHz/division  commencing  at  0.0  MHz. 
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Figure  4.7    Deconvolution  with  regularization;  d(k)  is  shown,  (b-f),  for  various 
input  and  output  signal  to  noise  ratios  and  degrees  of  filtering. 
SNR-x  75  dB,  SNR-y  46  dB :   (a)  h(k) ,  no  regularization,  n  ^  =   250  MHz, 
A  =  0;   (b)  moderate  regularization,  n  S7  =  100  MHz,  A  =  5;   (c)  strong 
regularization,  n  fi  =  74  MHz ,  A  =  6.   (d)  SNR-x  75  dB,  SNR-y  36  dB, 
strong  regularization,  n  SI  =  74  MHz,  A  =  6.   (e)  SNR-x  45  dB,  SNR-y 
46  dB,  moderate  regularization,  n  Q  =   100  MHz,  A  =  5.   (f)  SNR-x 
40  dB,  SNR-y  40  dB,  strong  regularization,  n  U  =  74  MHz,  A  =  6. 

Vertical  scale  2.5  x  10~-^/Div.  ,  Horizontal  scale  50  ns/Div. 
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Figure  4.8    An  example  of  the  imaginary  part  of  the  deconvolution  result. 
If  the  computation  procedures  of  the  DFT  (FFT)  were  exact, 
the  imaginary  part  would  vanish. 
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Figure  4.9    The  input  waveform  x(k)  and  the  magnitude  of  its 
DFT,  |x(n)|.   SNR-x   75  dB.   x(k)  is  shown  with 
an  abscissa  of  50  ns/div.  in  figure  5.15. 
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Figure   4.10  The   output   waveform,    y(k),    and   the  magnitude  of 

its   DFT,     |Y(n)|.      SNR-y   46   dB.      y(k)    is   shown 
with  an  abscissa   of   50  ns/div.    in   figure   5.15. 
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Figure  4.11  SNR-x  75   dB,    SNR-y   46   dB.      Deconvolution  without   regularization,    lH(n)|^g, 

h(k),    and   e(k).       |H(n)|  ,„    is    in   the   upper   right    hand   portion   of    the   upper 

dB 
figure;    vertical   scale   -10   dB/Div.    commencing  at   0.0   dB,    horizontal   scale 

25  MHz/Div.    commencing  at   0.0  MHz.      No   regularization  used;    this   relatively 
smooth  result    for  h(k)    is   possible   only  because  of   the   large  values   of 
SNR-x  and   SNR-y    (in  particular   SNR-y) .      Note   the   small  value  of   the   error 
function  e(k)    and   that   its   shape   is   comparable   to  y(k). 
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Figure   4.12  SNR-x   75   dB,    SNR-y  46   dB.      Deconvolution  with  moderate   regularization 


n  f^  =  100  MHz,  A  =  5, 
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|D(n)L„,  d(k),  and  e(k).   Note  the  filtering 
dB 


D(n)  j^  is  shown  in  the  upper  right  hand 
dB 


of  H(n)  to  obtain  D(n), 

portion  of  the  upper  figure,  Vertical  -  10  dB/Div.  commencing  at 

0.0  dB,  horizontal  25  MHz/Div.  commencing  at  0.0  MHz. 

50 


RE«.  PftRT  OF  T»C  OECONUOLUTION  RESULT 

Diy  ■  O.250E  -2 
mx  «  0  225E  -1 
MIN  — «.250E  -2 


HORIZONTAL  SOM^ 
OIU  «  288. ei  NS 


-1 

I 

/9=6 

\ 

> 

\, 

^ 

N 

\s 

s. 

•V 

\. 

\ 

-TrrJ 

n 

EXTENT  OF   "DC  USEFUL  SPECTRUM  - 

E!WR  FUtCTlCW 


DIU     «  0.227E  -1 

Tw:  «  0  112E   e 

PIN  —e.llSE  0 
^EPN  >  e  106E  -5 
SIOIO>  0.114E  -1 


74  mz 


■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

ff 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

m. 

.jii 

^^ 

Hi 

IM^U 

■i 

yyi 

m 

IH 

m 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

m 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

m 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

m 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

m 

EXTENT  OF  T«  USEFUL  SPECTRUM    -         74  WZ 


Figure   4.13  SNR-x   75   dB,    SNR-y   46   dB.      Deconvolution  with   strong  regularization 

n  ^  =   74   MHz,    A  =   6.       |D(n)|,    d(k),    and   e(k).      Note   the   filtering   of 
o 

H(n)    to   obtain  D(n).       |D(n)|         is   shown   in   the   upper   right   hand 

do 

portion  of  the  upper  figure,  vertical  -10  dB/Div.  commencing  at 
0.0  dB,  horizontal  25  MHz/Div.  commencing  at  0.0  MHz. 
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Figure   4.14  The   modified   output   waveform,    y(k),    and   the  magnitude   of   its 

a   SNR-y   of    36   dB 


DFT,     ,Y(n) I       ,    after   the   addition   of   10   dB   of   noise   to   yield 
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hand  portion  of  the  upper  figure,  vertical  scale  -10  dB/Dlv. 
commencing  at  0.0  dB,  horizontal  scale  25  MHz/Div. 
commencing  at  0.0  MHz. 
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Figure  4.16    SNR-x  75  dB,  SNR-y  36  dB.   Deconvolution  with  strong  regularization, 

n  ^  =  74  MHz,  A  =  6. 
o 

of  H(n)  to  obtain  D(n). 
portion  of  the  upper  figure.  Vertical  -10  dB/Div.  commencing  at  0.0  dB, 
horizontal  25  MHz/Div.  commencing  at  0.0  MHz. 
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Figure   4.17  The  modified   input   waveform,    x(k)    and   the  magnitude   of   its   DFT, 

|X(n)L    ,    after   the   addition  of   30   dB  of   noise   to   yield   a 
dB 

SNR-x  of   45   dB. 
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Figure   4.18  SNR-x,    45    dB,    SNR-y,    46   dB.      Deconvolution  without    regularization. 
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Figure   4.19  SNR-x   45   dB,    SNR-y   46   dB.      Deconvolution  with  moderate   regularization, 

n   fi  =    100  MHz,    A  =    5.       |D(n)ldB,    d(k),    and    e(k).      Note    filtering   of 


D(n)     ,„    is    shown    in   the   upper   right   hand   portion 
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H(n)  to  obtain  D(n)  , 

of  the  upper  figure,  vertical  scale  -  10  dB/div.  commencing  at  0.0  dB, 

horizontal  scale,  25  MHz/Div. ,  commencing  at  0.0  MHz. 
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Figure   4.21  The  modified   output   waveform,    y(k),    and   the  magnitude   of   its   DFT, 

|Y(n) I       ,    after   the   addition  of    6   dB   of   noise   to   yield   a   SNR-y   of 
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40  dB. 
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Figure  4.22  SNR-x  40  dB,    SNR-y  40  dB.      Deconvolution  without   regularization. 


|H(n)|^g,  h(k),  and  e(k). 


H(n)     is  in  the  upper  right  hand 
dB 


portion  of  the  upper  figure,  vertical  scale  -10  dB/Div.  conmencinj 
at  0.0  dB,  horizontal  scale  25  MHz/Div.  commencing  at  0.0  MHz. 
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Figure  4.23    SNR-x  40  dB,  SNR-y  40  dB.   Deconvolution  with  strong  regularization, 
n  fi  =  74  MHz,  A  =  6.   |D(n)L^,  d(k),  and  e(k).   Note  filtering 
of  H(n)  to  obtain  D(n).   |D(n)1^3  is  shown  in  the  upper  right  hand 
portion  of  the  upper  figure,  vertical  scale  -10  dB/Div.  coiranencing 
at  0.0  dB,  horizontal  scale  25  MHz/Div.  commencing  at  0.0  MHz. 
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Figure  4.24  The  time  domain  antenna  range  insertion  measurement  method. 
The  reference  input  waveform  x(k)  is  produced  by  the  system 
comprised  of  the  50  ohm  source  e(k) ,  the  source  antenna  a^ (k) , 
the  receiving  antenna  a_(k),  and  the  waveform  measurement 
system  s(k).   The  measurement  output  voltage  y(k)  is  produced 
when  a  (k)  replaces  a„(k). 
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Figure  4-31.   (a)  The  original  sequence,  f(k),  (dark  points)  is  N  units  in 
duration  and  is  constant  for  k  >  k  .   f (k)  is  extended 
by  N  units  (light  points). 

(b)  The  extended  f (k)  is  delayed  by  N  units  to  yield  f(k-N). 
f(k)  and  f(k-N)  both  consist  of  M  =  2N  points. 

(c)  f(k-N)  is  subtracted  from  f(k)  to  yield  f2(k). 
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Figure  4-32.   The  magnitude  of  the  spectra  for  F(n)  and  F  (m) 


The  discrete 


frequency  spacings  for  F(n)  and  F  (m)  are  fi  and  fi/2  due  to  the 
time  windows  NT  and  2NT,  respectively.   M=2N,  thus  the  total 
number  of  non-zero  values  in  F  (m)  is  the  same  as  in  F(n). 
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Figure  4.34   SNR-x  75  dB,  SNR-y  46  dB.   Deconvolution  with  moderate  regularization, 


Y 


100.   |D^(n)|,  d^Ck),  and  e(k).   Note  filtering  of  H  (n)  to  obtain 


D^(n).   |D^(n)|dB  is  shown  in  the  upper  right  hand  portion  of  the  upper 
figure,  vertical  -10  dB/Div.  commencing  at  0.0  dB,  horizontal  25  MHz/Div. 
commencing  at  0.0  MHz. 
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Figure  4.35      SNR-x   75   dB,    SNR-y  46   dB.      Deconvolution  with  strong  regular izatlon, 
Y  =   500.       [d    (n)l,    d    (k) ,    e(k).      Note   filtering  of  H^(n)    to  obtain 

|d    (n)|.       |d    (n)|dB    is    shown   in   the   upper    right   hand   portion   of    the 

upper    figure,    vertical   -10   dB/Div.    conmencing   at    0.0   dB,    horizontal 
25   MHz/Div.    commencing   at    0.0   MHz. 
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Figure   4.36.      SNR-x   75   dB,    SNR-y   46   dB.      Deconvolutlon  with  very   strong   regularization, 
Y=   10,000.       |D^(n)|,    d^Ck),    e(k).      Note   filtering   of   H^(n)    to   obtain  D    (n 
|D^(n)|dB   is    shown   in   the   upper   right   hand   portion  of    the   upper   figure, 
vertical,    -10   dB/div.    commencing   at   0.0   dB,    horizontal,    25  Mllz/div. 
commencing   at   0.0  MHz. 
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Figure  4.37   The  step-like  insertion  waveforms  x(k)  and  y(k),  and  their  duration 
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Figure  4.39  The  waveforms  and  spectral  magnitudes  for  very  strong 
regular izat ion,  y  =  10  . 
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5.   TIME  DOMAIN  DECONVOLUTION 

5.1  Introduction 

Time  domain  deconvolution  methods  operate  directly  on  the  time  domain  data  rather  than  on  transform 
domain  data.   There  are  certain  advantages  in  time  domain  operations.   Because  the  physical  manifestation 
of  the  signal  (pulse,  transient,  etc.)  occurs  with  the  evolution  of  time,  the  time  domain  is  the  natural 
or  physical  domain  of  the  signal.   Direct  operations  on  a  physical  signal  involve  operations  on  real 
quantities.   The  causal  and  stable  nature  of  a  signal  is  readily  recognized  in  the  time  domain  and  can 
be  directly  preserved  under  approximations  to,  and  operations  on  the  signal.   Also,  convolution  and 
deconvolution  processes  implemented  in  the  time  domain  only  employ  data  corresponding  to  times  such  that 
k  ;£  k  for  the  computation  of  y(k. )  and  h(k^),  respectively,  in 

y(k)   =  x(k)  *  h(k)  (5.1-1) 

h(k)   =  y(k)(l/*)  x(k)  (5.1-2) 

The  presence  of  noise  will  ultimately  yield  undesirable  results  (probably  unstable)  in  a  time  domain 
deconvolution  calculation;  this  fact  is  illustrated  later  in  the  discussion.   However,  filtering  or 
smoothing  of  the  time  domain  data  by  direct  operations  on  the  data  can  produce  stable  deconvolution 
results,  i.e.,  estimates  of  h(k).   The  apparent  promise  of  determining  h(k)  from  the  time  domain 
quantities  x(k)  and  y(k)  by  direct  operations  in  the  time  domain  is  a  sufficient  stimulus  to  try  and 
possibly  not  fail. 

The  time  domain  methods  presented  here  are  not  as  easy  to  apply  in  comparison  to  the  frequency 
domain  methods  presented  in  Chapter  4,  nor  are  the  results  as  smooth.   In  the  authors'  view,  the  time 
domain  methods  are  still  in  the  experimental  stages  of  development  while  the  frequency  domain  methods 
can  be  used  with  a  significantly  greater  degree  of  confidence  on  the  part  of  the  user.   However,  the 
promise  of  the  time  domain  method  remains  to  spur  further  development  studies. 

In  the  remaining  sections  of  this  chapter  the  following  topics  will  be  discussed: 

1.  Classical  time  domain  deconvolution 

2.  Iterative  time  domain  deconvolution 

3.  Regularization  or  filtering  in  the  time  domain 

4.  Experimental  results  of  classical  deconvolution  with  regularization  applied  to  the 
RG  58  C/U  transmission  line  problem  of  Chapter  4. 


5.2   Discrete  Classical  Time  Domain  Deconvolution 

5.2.1  Convolution 

The  discrete  classical  time  domain  deconvolution  method  is  an  algebraic  process  which  is  the  inverse 
of  the  convolution  algebraic  process  [26].   For  the  moment  consider  the  multiplication  of  the  two  poly- 
nomials f . (x)  and  f.(x). 

f,(u)   =  a  +  a,u  +  a„u^  +  •••  +  a^,  ^u^"-*-  (5.2.1-1) 

1         o    1     2  N-1 

f„(u)   =   b   +  b,u  +  b^u^  +  •••  +  b„  .u^"-*"  (5.2.1-2) 

I  O     1       /  N-1 


The  expression  for  the  coefficients  of  the  product 


f3(x)   =  f^(x)  f^Cx)  (5.2.1-3) 


is  identical  to  the  convolution  of  the  N  term  sequences 


h(k)   =   a^,  a^,  a^,  •••  ,  a^_^  (5.2.1-4) 

x(k)   =  b^,  b^,  b^,  •••  ,  h^_^  (5.2.1-5) 

and  is  given  by 

f 

y(k)  =  h(k)  *  x(k)   =  Coefficients  of  the  product,  f  (u)f„(u) 

k-1 
^   i?0  ""^^^"^  (5.2.1-6) 

In  the  notation  used  in  the  earlier  chapters,  the  convolution  result  for  y(k)  is  given  by 

k-1 
y(k)  =  ^    h(i)x(k-i)  (5.2.1-7) 

i=0 

Because  convolution  corresponds  to  multiplication  of  two  polynomials,  deconvolution  corresponds  to 
division  of  two  polynomials.   These  views  are  clearly  seen  in  the  Z-transform  representation  which  is 
discussed  in  section  5.2.2. 

For  the  causal  signals  h(k)  and  x(k),  both  are  equal  to  zero  for  k  <  0.   Expansion  of  (5.2-7)  for 
causal  signals  leads  to  a  result  in  which  the  k  =  0  term  is  zero;  hence  the  resultant  sequence  for  y(k) 
commences  with  k  =  1,  figure  5-la.   Accordingly,  deleting  the  initial  zero  values  h(0)  and  x(0), 
(5.2-7)  may  be  written  as 

k 
y(k)  =    X^hCi)  x(k+l-i)  (5.2.1-8) 

i=l 

which  yields  the  set  of  equations  for  the  y  sequence  values. 

yd)   =  h(l)  x(l)  (5.2.1-9) 

y(2)   =   h(l)  x(l)  +  h(2)  x(l)  (5.2.1-10) 

y(3)   =  h(l)  x(3)  +  h(2)  x(2)  +  h(3)  x(l)  (5.2.1-11) 

and  so  on.   Notice  that  y(k  )  is  computed  from  earlier  values  of  h(k)  and  x(k) ,  i.e.,  k  £  k  . 

5.2.2  Z-transform  Representation 

Some  of  the  properties  of  convolution  and  deconvolution  are  more  easily  developed  in  the  Z-transform 
domain  [29].   Before  proceeding  with  the  discussion,  some  of  the  properties  of  the  Z-transform  will  be 
briefly  described.   Given  a  pair  of  consistent  (exact)  discrete  waveforms  x(k)  and  y(k),  each  one  of  them 
can  be  transformed  to  the  complex  z  domain  by  the  Z-transf ormation 
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F(z)   =   Ef(k)  z  "  (5.2.2-1) 

k=0 


that  is,  f(k)  is  transformed  to  F(z);  hence 

CO 

X(z)   =   2]  x(k)  z"^  (5.2.2-2) 

k=0 

OO 

Y(z)   =   XI  y(k)  z'^  (5.2.2-3) 

k=0 

Furthermore,  the  discrete  convolution 

y(k)   =   x(k)  *  h(k) 
transforms  to 

Y(z)   =   X(z)  H(z)  (5.2.2-4) 

The  ratio  of  the  two  polynomials  for  Y(z)  and  X(z)  yields  H(z) , 


Z  y(k)  z  ^ 


„,  ,       k=0 

H(z)   = (5.2.2-5) 

E  x(k)  z~^ 
k=0 


which  is  another  polynomial  in  z, 


H(z)   =  ^      h(k)  z  ^  (5.2.2-6) 

k=0 

Consequently,  as  it  was  pointed  out  in  Section  5.2.1  the  coefficients  of  the  convolution  and 
deconvolution  sequences  are  generated  according  to  the  rule  for  polynomial  products  and  quotients, 
respectively. 


5.2.3  Deconvolution 

A  set  of  equations  for  generating  the  h  sequence  can  be  obtained  from  eqs.  (5.2.1-9)  through 
(5.2.1-11)  (and  so  on)  by  first  rearranging  them  to 

h(l)  x(l)  =  yd) 

h(2)  x(l)  =  y(2)  -  h(l)  x(l) 

h(3)  x(l)  =  y(3)  -  h(l)  x(3)  -  h(2)  x(2) 
and  so  on.   In  summation  form, 

k-1 

h(k)  x(l)   =   y(k)  -  22    h^i)  ^^^   +  1-i^  (5.2.3-1) 

i=l 

Dividing  x(k)  and  y(k)  by  x(l)  gives 

x(k)   =  jg^  (5.2.3-2) 
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y(k)   =  ^  (5.2.3-3) 

where 

x(l)  i     0  (5.2.3-4) 

x(l)   =  1  (5.2.3-5) 

Consequently,  h(k)  can  be  written  as 

_       k-1 

h(k)   =  y(k)  -   V  h(i)  x(k+l-i)  (5.2.3-6) 

1=1 

Briefly,  (5.2.3-6)  generates  the  following  sequence  for  h(k); 

h(l)   =  7(1)  (5.2.3-7) 

h(2)   =  7(2)  -  h(l)  x(2)  (5.2.3-8) 

h(3)   =  7(3)  -  h(l)  x(3)  -  h(2)  x(2)  (5.2.3-9) 

h(4)   =  7(^)  -  h(l)  x(4)  -  h(2)  x(3)  -  h(3)  x(2)  (5.2.3-10) 

h(5)   =  7(5)  -  h(l)  x(5)  -  h(2)  x(4)  -  h(3)  x(3)  -  h(4)  x(2) 

(5.2.3-11) 


Inspection  of  the  values  of  h(k)  on  a  line  to  line  basis  shows  that  the  earlier  values  of  the  impulse 
response  sequence,  e.g.,  h(l),  h(2),  ••',  h(k  -1),  are  used  to  calculate  the  specific  value  h(k  ). 
Because  the  h-values  can  be  expressed  in  terms  of  the  input  and  output  sequence  values,  x  and  y,  the  h 
values  are  actually  being  computed  from  the  x  and  y  values.   Note  that  h(l)  can  be  eliminated  from  the 
expression  for  h(2)  by  using  the  first  relation  for  h(l)  and  h(2). 

h(2)   =  7(2)  -  [7(1)]  x(2) 
Similarly,  for  h(3)  there  results, 

h(3)   =  7(3)  -  7(1)  x(3)  -  [7(2)  -  7(1)  x(2)]  x(2) 

continuing  the  process  obtains  the  partial  set  of  equations  for  h(k) , 

h(l)   =  7(1)  (5.2.3-12) 

h(2)   =   [7(2)  -  7(1)  x(2)]  (5.2.3-13) 

h(3)   =   [7(3)  -7(1)  x(3)]  -  [7(2)  -7(1)  x(2)]  x(12)  (5.2.3-14) 

h(4)   =   [7(4)  -  7(1)  x(4)]  -  [7(2)  -  7(1)  x(2)]  x(3)  -  [7(3)  -  7(1)  7(3)]  x(2) 

+  [y(2)  -  yd)  x(2)]  x^(2)  (5.2.3-15) 

h(5)   =   [7(5)  -  7(1)  x(5)]  -  [7(2)  -  7(1)  x(2)]  x(4) 

-  [7(3)  -  7(1)  x(3)[  x(3)  +  [7(2)  -  7(1)  x(2)]  x(2)  x(3) 

-  [7(4)  -  7(1)  x(4)]  x(2)  +  [7(2)  -  7(1)  x(2)]  x(3)  x(2) 

+  [7(3)  -  7(1)  x(3)]  x^(2)  -  [7(2)  -  7(1)  7(2)]  7^(2)  (5.2.3-16) 
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h(l)      = 

=      0 

h(2)      = 

=     a(l)    A(2) 

h(3)      = 

=     a(l)    A(3)    +  a(3)    A(2) 

and  so  on.   Defining  the  quantities, 

A(i)   =   y(i)  -  yd)  x(i)  (5.2.3-17) 

aCl)   =  1  (5.2.3-18) 

r-1      _ 
ct(r)   =  -   2  "(i)  x(r+l-i)  (5.2.3-19) 

i=l 

enable  h(k)  to  be  written  as 

(5.2.3-20) 

(5.2.3-21) 

(5.2.3-22) 

h(4)   =  a(l)  A(4)  +  a(2)  A(3)  +  a(3)  A(2)  (5.2.3-23) 

h(5)   =  a(l)  A(5)  +  a(2)  A(4)  +  a(3)  A(3)  +  a(4)  A(2)  (5.2.3-24) 

and  so  on.   Explicitly,  in  terms  of  (5.2.3-16)  through  (5.2.3-19),  h(k)  is  given  by 

k-1 
h(k)   =  yd)  +   2  a(i)  A(k+l-i)  (5.2.3-25) 

i=l 

k-1   •    _         _    _ 
=   yd)  +   22  "*^^^  [y(k+l-i)  -  yd)  x(k+l-i)]  (5.2.3-26) 

i=l 

which  is  the  general  result  for  h(k)  in  terms  of  x(k)  and  y(k). 

To  gain  further  insight  into  this  result  for  h(k),  (5.2.3-26),  consider  the  signals  shown  in 

figure  5.1.   x(k)  and  h(k)  are  given  signals  while  y(k)  is  calculated  from  x(k)  *  h(k),  figure  5.1a. 

Subsequently,  the  deconvolution  result,  h  (k) ,  is  calculated  from  y (k) (l/*)x(k) .   Note  that  both  y(k)  and 

h , (k)  commence  at  k  =  1  while  x(k)  and  h(k)  commence  at  k  =  0.   From  (5.2.3-17)  through  (5.2.3-19)  using 

the  data  of  figure  5.1,  there  results 

A(2)   =  7(2)  (5.2.3-27) 

(5.2.3-28) 
(5.2.2-39) 
(5.2.3-30) 

and  so  on.   Next, 

(5.2.3-31) 
(5.2.3-32) 

(5.2.3-33) 

-x(3)  +xh)]    x(2)  =  2x(2)  -  x-^(2)  (5.2.3-34) 

and  so  on.   Using  the  above  values,  calculate  the  deconvolution  of  y(k)  and  x(k)  to  obtain  h.(k).   To 

d 

insure  clarity,  the  deconvolution  result  for  h(k)  is  denoted  by  hj(k).   The  terms  in  h ,  (k)  will  be 

d  d 

uniquely  related  to  those  of  the  given  sequence,  h(k).   From  (5.2.3-20)  through  (5.2.3-24), 

k-1 
^d'^^^      "    E   a(i)A(k+l-i)  (5.2.3-35) 

i=l 

The  first  four  terms  are 
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A(3) 

= 

y(3) 

A(4) 

= 

7(4) 

A(5) 

= 

y(5) 

a(l) 

= 

1 

a(2) 

= 

x(2) 

a(3) 

= 

-x(3)    +  x2(2) 

a(4) 

= 

x(2)    x(3)    - 

h^(l)  = 

=      0 

h^(2)      = 

=     h(l) 

h^(3)      = 

=     h(2) 

h^(4)      = 

=     h(3) 

hjjd)   =  0  (5.2.3-36) 

h^(2)   =  a(l)A(2)  =  y(2)  (5.2.3-37) 

h^(3)   =  a(l)A(3)  +  cx(2)A(2)  =7(3)  -  x(2)  y(2)  (5.2.3-38) 

h^(4)   =  a(l)A(4)  +  a(2)A(3)  +  a(3)A(2) 

=  7(A)  -  x(2)  7(3)  +  [-x(3)  +  x^(2)]  7(2)  (5.2.3-39) 

Upon  substituting  the  values  of  y(k)  from  the  convolution  (5.2.1-8)  into  (5.2.3-36)  through  (5.2.3-39) 
and  recalling  that  x(k)  =  x(k)/x(l)  and  y(k)  =  y(k)/x(l),  the  following  values  for  h  (k)  are  obtained. 

(5.2.3-40) 
(5.2.3-41) 
(5.2.3-42) 
(5.2.3-43) 

Inspection  of  these  results  shows  that  h  (k)  is  equal  to  h(k-l),  figure  5.1b. 

Figure  5.2  shows  the  result  of  classical  deconvolution  using  two  noise  free  signals,  x(k)  and  y(k). 
In  this  case  x(k)  and  h(k)  were  specified  as 

x(k)  =  e  "^  (l-e~'^  );   k  =  0,1, •••,49  (5.2.3-44) 

h(k)  =  e'-^^^d-e'"'^) ;    k  =  0,1, •••,49  (5.2.3-45) 

which,  in  turn,  specified  y(k)  as 

y(k)  =  x(k)  *  h(k)     k  =  0,1, •••,49.  (5.2.3-46) 

The  data  for  y(k)  was  calculated  using  (5.2.1-8).   The  resultant  y(k)  was  used  with  x(k),  (5.2.3-44), 
in  (5.2.3-26)  to  calculate  the  deconvolution. 

h(k)  =  y(k)(l/*)  x(k)   k  =  0,1, •••,49  (5.2.3-47) 

In  applying  (5.2.3-6)  or  (5.2.3-26)  it  should  be  kept  in  mind  that  for  k  equal  to  zero,  h(k)  is  always 
zero  (causal  system);  consequently,  the  first  term,  zero,  must  be  added  to  the  sequence  for  h(k).   In 
the  present  example,  (5.2.3-6)  or  (5.2.3-26)  yields  a  sequence  of  49  terms,  and  including  zero  as  the 
first  term  in  the  sequence  gives  a  total  of  50  terms  in  the  sequence. 

5.2.4   The  Effect  of  Errors  in  Classical  Deconvolution 

The  effect  of  errors  or  noise  on  the  deconvolution  process  can  be  seen  by  considering  (5.2.3-6)  and 
its  partial  expansion  (5.2.3-7)  through  (5.2.3-11)  in  the  presence  of  errors.   For  the  moment,  recall 
the  deconvolution  example  given  above  in  Section  5.2.3.   The  noise  free  signals  x(k)  and  h(k) ,  (5.2.3-44) 
and  (5.2.3-45),  respectively,  were  used  to  calculate  a  noise  free  y(k),  (5.2.3-46);  this  y(k)  was  con- 
sistent with  x(k)  and  h(k)  through  the  convolution  equation  (5.2.1-8).   Now  imagine  that  one  value  of 
the  y(k)  sequence  is  changed  to  an  incorrect  value;  let  y(l)  be  altered  to  y(l), 

a. 

7(1)   =  7(1)  +  e  (1)  (5.2.4-1) 
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where  e  (1)  is  the  error  added  to  the  correct  yalue,  y(l).   Putting  y(l)  into  (5.2.2-7)  through  (5.2.2-11) 
yields  a  partial  set  of  the  erroneous  h(k)  s,  h(k), 

h(l)  =  h(l)  +  e  (1)  (5.2.4-2) 

h(2)  =  h(2)  +  e  (1)  [-x(2)]  (5.2.4-3) 

h(3)  =  h(3)  +e  (1)  [-x(3)  +x^(2)]  (5.2.4-4) 

h(4)  =  h(4)  +  e  (1)  [x(4)  -  2x(2)  x(3)  -  x-^(2)]  (5.2.4-5) 

h(5)  =  h(5)  +e  (1)  [x(5)  +x^(3)  +x^(2)  +x^(2)]  (5.2.4-6) 

and  so  on.   The  effect  of  the  error,  e  (1),  is  increased  with  increasing  k;  the  error  in  h(k)  is  simply 

e^(k)   =  h(k)  -  h(k)  (5.2.4-7) 

h 

Here,  the  set  of  errors,  e  (k) ,  in  h(k)  is  due  to  the  single  error  in  y(k) ,  e  (1).   In  practice,  the 
data  values  for  x(k)  and  y (k)  will  each  be  corrupted  by  noise;  consequently,  the  errors  in  h(k)  will  be 
dependent  upon  the  sets  of  error  for  both  x(k)  and  y(k).   Consequently,  one  car 
be  very  sensitive  to  the  errors  and  most  probably  will  be  an  unstable  sequence. 

The  analysis  can  be  extended  to  a  set  of  errors  in  be 
free.   Let  the  erroneous  values  of  y(k)  be  given  by  y(k). 


dependent  upon  the  sets  of  error  for  both  x(k)  and  y(k).   Consequently,  one  can  imagine  that  h(k)  will 
;ry  sensitive  to  the  errors  and  most  probably  will  be  an  unstable  sequence. 
The  analysis  can  be  extended  to  a  set  of  errors  in  both  y(k)  and  x(k) .   First,  let  x(k)  be  error 


y(k)   =   y(k)  +  ey(k)  (5.2.4-8) 

Putting  y(k)  into  the  general  result  for  h[x(k),  y(k)],  (5.2.3-26)  obtains 

h(k)   =  h(k)  +  e  (1)  +  y    :t(i)  [e  (k+l-i)  -  e  (1)  x(k+l-i),      (5.2.4-9) 

V  *.-  V  V 


y  i^i  y  y 


The  error  in  h(k)  is  e  (k)  and  is  given  by 


e^(k)   =   h(k)  -  h(k)  (5.2.4-10) 

h 

k-1  _ 

e,(l)  +  "Z,    "(i)  [e  (k+l-i)  -  e  (1)  x(k+l-i)]  (5.2.4-11) 


y 


i=l 


Next,  include  the  errors  in  x(k)  by  replacing  x(k)  and  a(k)  by 

x(k)   =  x(k)  +  e  (k)  (5.2.4-12) 


X 


% 


a(k)  =  a(k)  +  e^(k)  (5.2.4-13) 

The  final  result  for  the  erroneous  h  values  is  given  by 

_       k-1       _  _    _  k-1 

h(k)   =  yd)  +  J2    '-*(i)  [y(k+l-i)  -  yd)  x(k+l-i)]  -   2^  "(i)  yd)  e^(k+l-i) 
i=l  i=l 


k-1       _ 
y    e  (i){y(k+l-i)  -  yd)  [x(k+l-i)  +  e  (k+l-i)])  +  (con't,  next  page) 

i=  1  '' 
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k-1  _ 

+  e  (k)  +   X]  [a(i)  +  e  (l)]{e  (k+l-i)-e  (1) [x(k+l-i)  +  e  (k+l-i) ] }     (5.2.4-14) 


q    y 


The  error  in  h(k)  is  then 


e^(k)   =  h(k)  -  h(k)  (5.2.4-15) 

h 

k-1       _  k-1        _       _    _ 

e,(k)  -  Y  ^(i)  y(l)  e  (k+l-i)  +  Y    e„(i)  {y (k+l=i)-y (1) [x(k+l-i)+e  (1+1-i) ] } 


y 


i=l  i=l 


K-1  _ 

+  Yl     [■•(i)+e  (i)]  {e  (k+l-i)-e  (1)  [x(k+l-i)+e  (k+l-i)]}  (5.2.4-16) 

i=i  1  y  y  ^ 

To  illustrate  the  effect  of  noise  on  the  deconvolution  process,  noise  was  added  to  the  noise-free 
signals  (5.2.3-44)  and  (5.2.3-46)  so  that  each  signal  possessed  a  40  dB  signal  to  noise  ratio  (as 
defined  in  Section  4.1);  the  result  is  h(k) ,  (5.2.4-14),  and  is  shown  in  figure  5.3.   It  is  clear  that 
the  errors  (noise)  lead  to  an  unstable  result  which  bears  no  resemblance  to  h(k),  figure  5.2.   Similar 
results  are  also  obtained  for  the  iterative  deconvolution  method  which  will  be  discussed  in  section  5.3 


5.2.5  The  Critical  Dependence  of  Classical  Deconvolution 
on  the  First  Point  in  x(k). 

In  the  practical  application  of  the  classical  time  domain  method  to  the  deconvolution  problem,  the 
choice  of  the  first  term  in  the  input  sequence,  x(l),  is  critical.  If  x(l)  is  ill-chosen  the  resultant 
deconvolved  result  for  h(k)  will  not  converge.  To  gain  an  understanding  of  the  origin  of  this  feature, 
consider  the  Z-transf orms  of  the  two  known  sequences  x(m)  and  y(n)  having  their  durations  denoted  as 

M'  and  N' ,  respectively.   Thus, 

X(z)   =   x(l)  z"-"-  +  x(2)z~V  •••+  x(M') 

=  z~'''[x(l)  z^*  +  x(2)  z^'"-^+  •••+  x(M+l)];   M'-1  =  M  (5.2.5-1) 

Y(z)   =   yd)  z"-^  +  y(2)  z"^  +  •••  +  y(N') 

=   z~^[y(l)  z^  +  y(2)  z^""""  +  •••  +  y(N+l)]   N'-l  =  N  (5.2.5-2) 


and 


respectively,  where 


N>M  (5.2.5-3) 

The  Z-transf orm  of  the  impulse  response  corresponding  to  x(k)  and  y(k)  is 

nil)      =     Y(n)/X(m)  (5.2.5-4) 

When  N  equals  M,  H("".)  becomes  H(n)  and  is  a  constant. 

H(n)   =  ^^  =  constant  =C;   N  =  M  (5.2.5.5) 

A(,n) 

The  inverse  Z-transf orm  of  H(n)  is  then  the  impulse  sequence  C,  0,  0,  •••  0  in  which  the  first  term  is  C 
and  is  followed  by  N-zeroes. 

By  synthetic  division  H(z)  is  given  by 
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H(z)  =  -^!^z^-.ry(2)-^^^%^1^ 


N-M-1 


N-M  x(l)        \  y^^'  5^(2^)   J   x(l) 

z     '  L 


i^^n^  _  y(i)x(3)1  _  r  .   _  y(i)x(2)1  2L(2i  j  z^-^-^   ..  j 

j[y^^^     x(l)   J    [y(2)      ^^^^  J  ^(^)  j    ^(^)  +    ^ 


(5.2.5-6) 
Using  the  variables  x(k)  and  y(k)  which  are  the  x(k)  and  y(k)  normalized  to  x(l),  H(z)  can  be  written 
as 

H(z)  =  7(1)  +  [7(2)-7^(l)x(2)]z"^  +  {7(3)-7(l)x(3)-[7(2)-7(l)x(2)]x(2)}z"^  +  • 

(5.2.5-7) 
Now  consider  that  the  point  x(l)  is  corrupted  by  an  error,  e  ,  and  thus  changed  to  x(l) , 

x(l)   =  x(l)  +  e  (5.2.5-8) 

while  all  other  values  of  x(k)  and  y(k)  are  correct.   The  normalized  x(m)  and  y(n)  would  be  given  by 

-/  ^      1     x(2)      _.      x(M') 

x(.m)   =   1,  — 7TT7 >     »  — 7T^C~Z 

x(l)+  e         x(l)  +  e 

X  X 

=   1,  x(2)  +  e  (2),  •••,  x(M')  +  e  (M')  (5.2.5-9) 

-(   ^  =    y(l)         y(2)    ...    y(N') 

^^    '  x(l)  +  e    '   x(l)  +  e   '     '  x(l)  +  e 

XX  X 

=   yd)  +  e  (1),  y(2)  +  e  (2),  •••,   y(N')  +  e  (N')  (5.2.5-10) 

respectively.   Then  the  only  correct  term  in  x(m)  would  be  x(l)  which  equals  x(l)  =  1.   The  y(n) 
sequence  would  be  scaled  to  the  erroneous  value  x(l)  +  e  .   Consequently,  the  ratio  for  H(z)  would  also 
be  erroneous, 

J^(z)  =  li^  =   7(1)  +  [7(2)  -  7(1)  x(2)]  z  "  +  {  y(3)  -  yd)  x(3) 

X(z) 

-      [7(2)  -  yd)  x(2)]  x(2)}  z    +•••  (5.2.5-11) 

The  inverse  Z-transform  of  H(z)  gives  the  erroneous  sequence  h(k),  (5.2.4-14),  and  an  error  in  h(k)  of 
,  (5.2.4-16). 
Consequently,  it  has  been  shown  that  the  error  terms,  i.e., 

e^(k)  :  e^(2),  e^(3) ,  ••• 
-a^^)  =  ^a(2),  e^(3),  ••• 
e  (k)  :   e  (1),  e  (2),  ••• 

y      y    y 
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e^(k),  (5.2.4-16) 
h 


e^(k)  :   e^(2),  e^(3),  •.• 

have  been  generated  by  the  single  error  in  x(l)  through  the  normalization  processes  (5.2.5-9)  and 
(5.2.5-10).   Thus  the  selection  of  x(l)  is  critical  to  the  deconvolution  process. 


5.3   Discrete  Iterative  Time  Domain  Deconvolution 

5.3.1   Iterative  Method 

By  using  an  iteration  process  it  is  also  possible  to  solve  the  deconvolution  equation  by  an 
iterative  convolution  process  [4,5,27,28]  proposed  by  P.  H.  van  Cittert.  [28].   The  principle  of  this 
method  is  to  form  successive  approximations,  »h.(k),  to  the  unknown  system  impulse  response,  h(k) ,  using 
the  convolution  equation, 

y(k)   =  h(k)  *  x(k)  (5.3.1-1) 

For  physically  realizable  signals  the  duration  of  x(k)  is  less  than  y(k).   As  a  first  approximation 
to  h(k),  h^(k),  use  y(k), 

h^(k)   =  y(k)  (5.3.1-2) 

Using  h^ (k)  in  the  convolution  equation  (5.3.1-1)  gives  the  result  y  (k) . 

y^(k)   =  h^(k)  *  x(k)  (5.3.1-3) 

y.  (k)  differs  from  the  true  result  y(k)  by  the  error,  [Ay(k)] 

[Ay(k)]^   =   y(k)  -  y^(k)  (5.3.1-4) 


which  corresponds  to  an  error  in  h(k)  of 


[Ah(k)]^  =   [Ay(k)]^  (1/*)  x(k) 


(5.3.1-5) 


correcting  h  (k)  by  adding  [Ah(k)]   yields  a  second  approximation  to  h(k),  h„(k). 

h2(k)   =  h^(k)  +  [Ah(k)]^  (5.3.1-6) 

=  h^(k)  +  [y(k)  -  h^(k)  *  x(k)]  (5.3.1-7) 

Consequently,  the  i-th  approximation  is  related  to  (i-l)th  approximation  by 

h.(k)   =  h^_^(k)  +  [y(k)  -  h^_^(k)  *  x(k)]  (5.3.1-8) 

When         ^ 

h.(k)  *   x(k)   =  y(k),  (5.3.1-9) 


1 
then 


h^(k)   =  h^^^(k)   =  h(k)  (5.3.1-10) 


Here,  the  y(k)  and  x(k)  are  considered  to  be  noise  free.   Of  course,  they  are  not  in  practice,  and  the 
noise  does  limit  the  direct  application  of  the  iterative  method. 

Continuing  the  discussion  on  the  iterative  method,  apply  the  Z-transformation  to  (5.3.1-8). 
The  result  is 

Y(z)  -  H^(z)  X(3)  +  H^(z)   =  H._|_^(z)  (5.3.1-10) 

The  form  of  (5.3.2-7)  parallels  that  of  (5.3.1-8)  with  the  convolution  operation  replaced  by  a  product. 
The  iteration  process  starts  with  the  first  estimate  for  H(z)  being  chosen  as  Y(z), 

H^(z)   =   Y(z)  (5.3.1-11) 

Putting  H  (z)  into  (5.3.2-4)  yields  for  U^(z) 

n^(z)      =  Y(z)[2-X(z)]  (5.3.1-12) 

Repeating  the  process  for  the  next  values  yields 


n^iz)      =   Y(z)  +  H2(z)[l  -  X(z)]   =   Y(z)  {1  +  [2  -  X(z)][l  -  x(z)]} 


(5.3.1-13) 


H,(z)   =   Y(z)  {1  +  [l-X(z)]  +  [2-X(z)][l-Y(z)]^}  (5.3.1-14) 

4 

H^(z)   =  Y(z)  {1  +  [l-X(z)]  +  [l-X(z)]^  +  [2-X(z)][l-X(z)]^}        (5.3.1-15) 

H,(z)   =   Y(z)  {1  +  [l-X(z)]  +  [l-X(z)]^  +  [l-X(z)]^ 
6 

+  [2-X(z)][l-X(z)]^}  (5.3.1-16) 


and  so  on.   The  result  for  H.(z)  is  a  power  series  of  the  form 


[±  +  B^   X(z)    +  B2   X(z)^   +    •••    +  B._^(z)'-    -^J 


H^(z)      =      Y(z)   1  i  +  B,    X(z)    +  B^    X(z)      +    •••    +  B._, (z)  \  (5.3.1-17) 

whose    inverse   Z-transform   is 

h.(k)      =      iy(k)    +  B^y(k)    *   X(k)    +  B2y(k)    *  X(k)    *  X(k)   +   ••• 

+  B.    ,    y(k)    '•=  x(k)    *  x(k)    ■"    •••x(k)  (5.3.1-18) 

1-1   " 

Hence,  it  is  seen  that  the  direct  discrete  iterative  method  requires  a  large  number  of  convolution 

operations  which  in  turn  will  require  a  large  amount  of  computation  time. 

Using  the  x(k)  and  y(k)  of  (5.2.3-44,  5.2.3-46),  the  direct  discrete  iterative  deconvolution  method 
was  employed  to  calculate  h(k);  the  results  for  500  and  1000  iterations  are  shown  in  figures  5.4  and  5.5 
where  it  is  seen  that  the  latter  result  is  identical  to  the  result  obtained  by  the  direct  discrete 
classical  deconvolution  method.   However,  the  discrete  iterative  method  required  45  minutes  of  computa- 
tion time  as  compared  to  10  seconds  for  the  classical  method. 


5.3.2  The  Effect  of  Errors  in  Iterative  Deconvolution 

If  y(k)  contains  an  error,  say  for  k  =  0,  then  the  erroneous  output  sequence  y(k)  is 

y(k)   =  6^(0)  +y(0),  yd),  y(2),  •••  (5.3.2-1) 

where  e  (o)  is  the  error  in  y(0).   The  Z-transform  of  y(k)  is 

Y(z)   =  e  (0)  +  y(0)  +  yCDz""""  +  y{2)  z"^+  •••  (5.3.2-2) 

=   e  (0)  +  J2      y(k)  z"^  (5.3.2-3) 

■^       k=0 

=  e  (0)  +  Y(z)  (5.3.2-4) 

Using  the  erroneous  output  Y(z)  in  (5.3.1-17)  yields  the  erroneous  i-th  approximation  to  the  deconvolution 
result 

H^(z)   =   [e  (0)  +  Y(z)]  [i  +  B^X(z)  +  B^X^  {z)    +    •••  +  B  ._^x'^"-'-(z)  ]   (5.3.2-5) 

=  +  H.(z)  +e  (0)  F.[X(z)]  (5.3.2-6) 

where  the  function  of  X(z)  is  defined  by 

2  ^-1 

F^[X(z)]   =   [i  +  B^X(z)  +  B2X  (z)  +  •••  +  B^_^X(z)]  (5.3.2-7) 

Consequently,  the  effect  of  the  error  e  (0)  increases  with  increasing  k.   This  result  is  similar  to  the 
error  result  for  the  classical  deconvolution  technique  with  a  single  error  (5.2.4-6). 

For  a  set  of  errors,  e  (k) ,  in  y(k),  i.e.,  for  an  erroneous  output  sequence  of  the  form 

y(k)   =  e  (k)  +  y(k)  (5.3.2-8) 

the  Z-transform  of  the  sequence  is 

Y(z)   =  E  (z)  +  Y(z)  (5.3.2-9) 

Using  (5.3.2-9)  in  (5.3.1-17)  yields  the  erroneous  i-th  approximation  to  the  deconvolution  result 

H^(z)   =   [E  (z)  +  Y(z)]  F^[X(z)]  (5.3.2-10) 

=  H.(z)  +  E  (z)  F.[X(z)]  (5.3.2-11) 

whose  inverse  Z-transform  yields  the  erroneous  time  domain  sequence,  h.(k), 

h.(k)   =  h.(k)  +  e  (k)  *  f.(k)  (5.3.2-12) 

where  f . (k)  is  the  inverse  Z-transform  of  F.[X(z)].  h.(k)  is  the  y-noise  or  error  corrupted  result  for  the 
noise  free  i-th  iteration  approximation,  h.(k),  to  the  system  impulse  response,  h(k). 

If  x-noise  or  errors,  e  (k) ,  are  also  present  they  transform  to  E  (z) ;  and  then  X(z)  becomes  X(z) , 
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X(z)   =  X(z)  +  E^(z)  (5.3.2-13) 

and  the  corrupted  i-th  iteration, is  given  by 

H.(z)   =  H.(z)  +  Ey(z)  G.[X(z)  +  E^(z)]  (5.3.2-14) 

G.[X(z)  +  E  (z)]  is  shown  in  expanded  form  in  (5.3.2-16).   The  error  in  H.(z)  is  then 
x^         X  X 

Ej^^(z)   =   H^(z)  -  H.(z) 

=  E  (z)G  [X(z)+E  (z)]  (5.3.2-15) 

y    1       X 

=   E  (z)  ]  i  +  B  [X(z)+E  (z)]  +  B„[X(z)+E^(z)]V  •••+B.  , [X(z)+E  (z)]^"^} 
y    '      -L       X        z       X  i—x       X 

(5.3.2-16) 

In  the  time  domain,  the  i-th  corrupted  iteration,  h.(k),  for  h(k)  is  given  by  the  inverse  Z-transform 

of  (5.3.2-14). 

h^(k)   =  h^(k)  +  Ey(k)  *  g^[x(k),  e^(k)]  (5.3.2-17) 

where  g.(k)  is  the  inverse  Z-transform  of  G.[X(z)  +  E  (z)]. 

Figure  5.5  shows  the  result  for  h.(k)  with  i=750  (iterations)  using  noise  corrupted  values  of 

•^  """% 

(5.2.3-44)  and  (5.2.3-46)  for  x(k)  and  y(k),  respectively, 

x(k)   =  x(k)  +  e  (k)  (5.3.2-18) 

y(k)   =  y(k)  +  e  (k)  (5.3.2-19) 

where  e  (k)  and  e  (k)  were  such  that  the  SNR  of  x(k)  and  SNR  of  y(k)  were  both  40  dB.   Note  that  the 
x         y 

result  is  unstable  and  bears  no  resemblance  to  the  system  response  h(k),  (5.2.3-45).   Figure  5.5  is  very 
similar  to  figure  5.2,  the  result  for  noise  corrupted  classical  deconvolution  under  the  same  x  and  y 
signal  to  noise  ratios  (40  dB) . 

5.4   Time  Domain  Regularization 

5.4.1  Averaging  Filter 

In  sections  5.2  and  5.3  two  direct  methods  of  deconvolution  were  presented.   By  direct  it  is  meant 
that  the  mathematical  method  is  directly  applied  to  the  data  x(k)  and  y(k).   As  we  demonstrated  in 
sections  5.2  and  5.3  the  presence  of  noise  or  errors  in  the  data  destroys  the  consistency  between  x(k) 
and  y(k);  that  is  to  say,  the  computed  data  x(k)  and  y(k)  are  no  longer  mathematically  related  to  each 
other  through  the  impulse  response  of  the  system,  h(k) .   Thus,  the  direct  application  of  the  mathematical 
methods  to  the  computed  data  x(k)  and  y(k)  leads  to  an  erroneous  and  divergent  or  unstable  result,  h(k) . 
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Both  the  classical  and  iterative  methods  applied  to  noisy  signals  lead  to  the  same  type  of  Incorrect 
and  divergent  results. 

By  filtering  the  noisy  waveforms  x(k)  and  y(k)  both  in  the  same  prescribed  way,  the  functions  x(k) 
and  y(k)  can  be  smoothed  or  regularized  so  that  when  the  mathematical  deconvolution  method  is  applied  to 
them,  the  process  yields  a  stable  convergent  result,  d(k),  which  approximates  h(k).   The  smoothing  is 
accomplished  by  a  regularization  filter.   Figures  5.7  and  5.8  show  the  resultant  estimates  for  h(k) 
obtained  by  first  filtering  the  corrupted  signals  x(k)  and  y(k)  which  were  used  in  earlier  examples. 
Figure  5.7  was  obtained  by  using  the  filtered  signals  in  the  discrete  classical  deconvolution  method, 
(5.2.3-6)  or  (5.2.3-26).   Simarily,  figure  5.8  was  produced  by  the  regularized  discrete  iterative  decon- 
volution  method,  i.e.,  filtered  x(k)  and  y(k)  in  equation  (5.3.1-8).   Comparison  of  figures  5.7  and  5.8 
shows  that  the  two  methods  yield  very  similar  results,  but  again,  the  computation  time  for  the  iterative 
method  was  45  minutes  compared  to  10  seconds  , for  the  classical  deconvolution  method.   The  notation  d(k) 
denotes  an  estimate  of  the  true  impulse  response  h(k)  which  was  estimated  from  the  filtered  or  regularized 

data  x(k)  and  y(k). 

^\,  'X, 

The  regularizzation  filter  applied  to  x(k)  and  y(k)  was  a  simple  one,  a  four-point  averaging  filter: 

the  value  of,  say  x(k) ,  at  the  point  k  was  set  equal  to  the  average  or  mean  value  of  the  values  of 

x(k) ,  x(k+l),  x(k+2),  and  x(k+3);  i.e.,  x(k)  was  set  equal  to  x(k) . 

N 

^(k)   =  ^   E  ^(J)  (5.4.1-1) 

j=k 

where  N  =  k  +  3.  Without  the  filter  (5.4.1-1),  the  S/N  ratios  for  x  and  y  of  40  dB  led  to  failure  of 
the  direct  deconvolution  methods.  With  the  filter  (5.4.1-1),  the  S/N  ratios  of  20  dB  led  to  failure; 
consequently,  it  could  be  said  that  there  was  a  20  dB  enhancement  of  the  deconvolution  process. 

A  major  task  in  time  domain  deconvolution  is  the  synthesis  of  a  suitable  regularization  filter. 
The  example  described  above  employed  the  simple  averaging  filter  (5.4.1-1);  however,  to  achieve  improved 
estimates  of  the  impulse  response,  d(k),  more  sophisticated  regularization  filters  need  to  be  synthesized. 
In  the  next  section  a  more  sophisticated  and  flexible  filter  is  described. 


5.4.2   Composite  Filtering 

5.4.2.1   Introduction 

In  the  last  section  an  averaging  filter  was  briefly  discussed  and  applied  to  two  deconvolution 
problems.   The  averaging  filter  is  equivalent  to  a  system  having  a  rectangular  impulse  response,  h  (k) . 


h  (k) 


1'   0  1k  <k^  (5.4.2.1-1) 


0,   k  <  k  <  N-1 


In  other  words,  the  impulse  response  h  (k)  consists  of  N  discrete  values,  k^  +  1  of  which  are  equal  to 

a  1 

unity,  while  the  remaining  values  are  zero. 

If  a  given  waveform  sequence  x(k)  is  convolved  with  h  (k)  using  (5.2.1-7)  or  (5.2.1-8),  the  result 
is  the  same  as  that  produced  by  the  averaging  operation,  (5.4.1-1). 

Presented  below  is  a  filter  consisting  of  three  operations  which  are  (1)  averaging,  (2)  least-square 
first-order  polynomial  fitting,  and  (3)  Gaussian  weighting  or  smoothing.  Each  one  of  these  operations  is 
represented  by  an  adjustable  parameter.   Consequently,  the  overall  composite  filter  possesses  three 
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parameters,  L,  S,  and  A  corresponding  to  averaging,  fitting,  and  weighting,  respectively. 

5.4.2.2   Filter  Description 

The  main  characteristic  of  this  filter  is  that  it  is  possible  to  adjust  its  transfer  characteristic 
over  a  relatively  wide  range  by  adjusting  the  three  parameters  L,  S,  and  A.   The  transfer  properties  will 
be  considered  in  the  next  section  while  here  the  filter  structure  will  be  presented.   L  and  S  are  integers 
while  A  is  a  positive  real  number.   Let  W(k)  be  the  discrete  waveform  to  be  filtered. 

L  is  the  running  averaging  parameter  and  is  equal  to  the  time  window  (or  number  of  points)  of  the 
averaging  operation.   Each  value  of  W(k) ,  say  W(k  ),  is  replaced  by  the  mean  value  of  W(k. )  and  the  L-1 
values  of  W(k)  following  W(k  ) ,  i.e., 


New  W(k)   =  L   S  W(k+i)  (5.4.2.2-1) 


S  is  the  polynomial  fitting  parameter  and  is  equal  to  the  time  window  (or  number  of  points)  over 
which  the  least  square,  first  order  polynomial  fitting  [30]  is  performed.   The  window  is  shifted  along 
the  waveform  in  such  a  way  that  each  point  of  the  waveform  W(k) ,  say  W(k. ) ,  comes  S  times  in  the  time 
window.   For  each  position  of  the  window  the  S  selected  points  are  replaced  by  a  new  set  of  S  values 
deduced  from  the  fitting.   The  first  time  the  point  W(k  )  appears  in  the  window  it  is  transformed  to  a 
new  value  W(k  ,1);  the  second  time  the  value  W(k  ,2)  is  obtained  and  so  on  until  the  value  W(k  ,S)  is 
obtained. 

A  is  the  characteristic  constant  of  a  Gaussian  weighting  function  which  is  applied  to  the  S  values 
of  W(k  ,  i) , 

S 
New  W(k)   S  W  ^(k)   =   ^   a(i)  W(k,i)  (5.4.2.2-2) 

'         i=l 


where 


i=l 
The  expanded  form  of  (5.4.2.2-2)  is  then 


_.2,  2 
a(i)   =   a  e  "■  ^  (5.4.2.2-3) 

S 
^   a(i)   =  1  (5.4.2.2-4) 


-2  -2  2  -2 

Wg  ^(k)   =   a  I  e"'^   W(k,l)  +  e""*^  W(k,  2)  +  •  •  •  +  e"^  "^  W(k,S)|     (5.4.2.2-5) 


Te  ^        W( 
=  le     +e      +'"+e 


where 

-2      ,.-2  o2.-2 


(5.4.2.2-6) 


Summarizing,  the  composite  filter  employs  three  filtering  or  smoothing  techniques:  averaging,  least 
squares-first  order  fitting,  and  Gaussian  weighting,  with  the  corresponding  adjustable  filter  parameters 
L,  S,  and  A,  respectively.  The  filter  is  divided  into  two  filters,  the  L-averaging  filter  and  the  S,  A- 
fitting/weighting  filter. 
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5. A. 2. 3  Frequency  Domain  Analysis  of  the  Composite  Filter 

Frequency  domain  analysis  of  the  composite  time  domain  filter  (or  any  other  time  domain  filter)  is 
not  necessary  for  the  application  of  time  domain  filtering.  It  is  presented  here  to  provide  the  reader 
with  additional  insight  to  the  filtering  processes  by  establishing  a  connection  to  the  frequency  domain 
and  thus  to  some  conventional  frequency  domain  concepts  and  terminology. 

Upon  applying  each  filter,  (5.4.2.2-1)  and  (5.4.2.2-2),  to  an  impulse  (a  single  data  point)  the 
impulse  response  of  each  filter  can  be  obtained.   The  resultant  impulse  response  sequences  can  then  be 
Fourier  transformed  to  the  frequency  domain. 

The  transfer  function  of  the  averaging  (L-parameter)  filter  (5.4.2.2-1)  is  characterized  by  the  cut- 
off frequency  fco  whose  value  decreases  as  L  increases.   For  example,  for  N  =  1024  points  and  a  total 
time  window  of  200  ns  the  corresponding  values  for  L  and  fco  are  as  follows:   L=  2,  3,  4,  5,  6,  7,  8,  9 
and  10,  and  fco  (MHz)  =  500,  340,  255,  205,  170,  145,  130,  115  and  105,  respectively. 

The  transfer  function  for  an  L-f liter  with  L  =  3  is  shown  in  figure  5.9A.   Also  shown  is  the 
transfer  function  for  an  S,  A  filter  with  S  =  5,  A  =  1  and  the  cascade  (logarithmic  sum)  of  the  two 
filters.   Note  that  the  L  filter  cut-off  frequency  was  chosen  to  occur  at  the  first  maximum  of  the  S,  A 
filter  to  yield  a  sharper  cut-off  characteristic  and  increased  attenuation  above  fco. 

The  transfer  properties  of  the  S,  A  filter  depends  upon  the  two  parameters  S  and  A.   Figure  5.10 
shows  how  the  value  of  A  affects  the  transfer  properties  with  S  being  held  constant  and  equal  to  5. 
Similar  dependencies  on  A  appear  in  figures  5.11  and  5.12;  however,  figure  5.10  provides  a  systematic 
view  of  the  effects  of  A. 

In  figure  5.10,  A  is  varied  over  the  range  of  0.5  through  5.0.   For  S  =  5,  the  value  of  A  =  2.05 
is  of  special  significance.   As  A  is  increased  from  0.50,  the  first  two  minima  shift  towards  each  other 
and  merge  together  when  A  equals  2.05.   The  resultant  cut-off  characteristic  or  minimum  is  very  deep, 
about  100  dB.   As  A  is  increased  above  2.05,  the  cut-off  frequency  remains  the  same  while  the  overall 
transfer  function  broadens  out  and  becomes  smoother. 

In  figures  5.11  and  5.12  examples  of  composite  filtering  are  shown.   Figure  5.13  illustrates  the 
S,  A  filter  cut-off  frequency  dependence  on  the  values  of  S  and  A,  and  can  be  used  as  a  design  tool  for 
selecting  S  and  A  to  achieve  a  given  cut-off  frequency.   The  cut-off  frequency,  fco,  can  be  scaled 
with  Af;  in  the  present  case  Af  =  50  MHz/div. 

In  figure  5.13  for  a  given  fco  there  are  multiple  choices  for  the  pair  of  S,  A  values.   For  example, 
if  fco  is  desired  to  be  120  MHz,  figure  5.13  indicates  that  two  possible  choices  for  S,  A  are  12,  4.4 
and  13,  5.61.   Some  contours  of  limited  extent  are  shown  for  A  equal  to  1.0  through  6.0  along  with 
single  points  for  A  over  the  range  of  2.51  through  6.49. 

For  1.0  ^  A  2l  2.5,  the  cut-off  frequencies  of  the  S,  A  filters  in  figures  5.9  and  5.10  can  be 
compared  to  those  in  figure  5.13.   Also  by  using  a  time  window  of  200  ns  (Af  =  50  MHz)  in  figures  5.11 
and  5.12,  the  cut-off  frequencies  in  figures  5.11  and  5.12  can  be  compared  to  those  in  figure  5.13.   For 
example,  from  figure  5.11  fco  =  250  MHz  occurs  for  S,  A  equal  to  6,  2.40.   Inspection  of  figure  5.13 
shows  that  for  S  =  6,  A  would  be  between  2  and  3. 
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5. A. 3   Application  of  Time  Domain  Filtering 
5.4.3.1   Practical  Procedure 

Time  domain  deconvolution  is  a  progressive  process  (one  point  after  another) .   Knowledge  or  selection 
of  the  first  point  in  the  input  sequence,  x(l),  is  critical  as  was  shown  in  section  5.2.5.   The  selection 
o'f  x(l)  is  accomplished  by  trial  and  error  using  the  deconvolution  algorithm;  however,  this  is  not  as 
random  an  exercise  as  it  may  appear  to  be.   The  search  for  x(l)  is  performed  using  only  a  small  part  of 
the  available  data  points,  say  50  points  out  of  1024.   Typically,  a  choice  is  made  for  x(l)  and  the 
deconvolution  algorithm  is  run  for  50  points  or  so.   If  the  result  is  stable,  then  the  algorithm  is 
applied  to  increasingly  larger  numbers  of  data  points  until  the  entire  data  span  is  Included,  with  x(l) 
being  adjusted  until  a  value  is  found  which  provides  an  overall  stable  deconvolution  result. 

There  are  two  factors  that  strongly  interact  with  the  selection  of  x(l) :   (1)  data  offset  error  and 
(2)  noise  in  the  x  and  y  data.   The  offset  error  correction  can  be  applied  before  commencing  the  decon- 
volution operation  (refer  to  section  4.6.3).   Noise  in  the  data  can  only  be  eliminated  by  filtering; 
consequently,  some  combination  of  offset  correction,  x(l)  selection,  and  filtering  is  necessary  to 
achieve  an  acceptalbe  deconvolution  result.   All  of  this  is  accomplished  in  an  iteractive  mode  of 
operation.   First,  the  operator  selects  the  offset  correction,  x(l),  and  the  filter  applied  to  x(k)  and 
y(k),  and  then  judges  the  quality  of  the  displayed  deconvolution  result.   The  step  by  step  procedure  is 
as  follows: 

1.  Preparation  of  the  waveforms,  x(k)  and  y(k). 

a.  Offset  correction 

b.  Choose  x(l) 

2.  Filtering 

a.  Choose  the  three  filter  parameters,  L,  S  and  A 

b.  Apply  the  filter  to  both  waveforms,  x(k)  and  y(k) 

3.  Deconvolution 

a.  Is  the  result  stable,  i.e.,  does  it  converge? 
If  "no"  go  back  to  step  1  or  step  2. 

b.  Is  the  result  smooth? 

If  "no"  go  back  to  step  2 

c.  Continue 

4.  Error  Evaluation 

a.  Is  the  error  acceptable? 
If  "no"  go  back  to  step  2 

b.  Continue 

5.  Print  out  results:   d(k)  with  related  quality  parameters  (mean  error,  standard  deviation,  etc.). 
d(k)  is  the  estimate  of  the  impulse  response. 

The  development  of  the  appropriate  filter  (in  2a  above)  should  nominally  adhere  to  the  following 
rules : 

1.   Start  with  the  filter  L  =  2 ,  S  =  3,  A  =  100 
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2.  Keeping  L  =  2,  and  A  =  10,  Increase  S  until  the  result  is  stable  (converges). 

3.  Keeping  S  and  A  constant  increase  L  until  the  result  is  smooth. 

4.  After  L  and  S  are  determined,  compute  the  error  criteria  (section  3.4). 

-  if  the  mean  value  of  the  error,  e,  is  greater  than  zero,  decrease  the  value  of  A. 

-  if  the  mean  value  of  e  is  less  than  zero,  increase  S  until  e  becomes  greater  than  zero. 
The  smaller  the  value  of  e,  the  better  is  the  result. 

5.  When  e  is  close  to  zero,  compute  the  standard  deviation,  0. 

-  if  o  is  too  large,  increase  L  or  S  and  adjust  A  to  reduce  O. 

In  summary,  as  it  is  evident  from  the  above  described  procedures,  time  domain  deconvolution  is  a 
progressive  process.   The  programs  used  to  implement  time  domain  deconvolution  were  written  to  work  in 
an  iteractive  mode.   Knowing  that  the  result  must  be  convergent,  smooth,  and  close  to  the  true  solution 
(small  error  function),  the  operator  selects  the  offset  correction,  x(l) ,  and  the  filter,  judges  the 
quality  of  the  displayed  result,  and  then  modifies  those  selections  until  an  acceptable  result  is 
obtained. 

5.4.3.2  Estimation  of  the  Insertion  Impulse  Response 
of  300  Meters  of  RG  58  C/U  Coaxial  Cable 

The  following  examples  of  time  domain  deconvolution  employs  the  same  x(k)  and  y(k)  data  as  were  used 
in  the  frequency  domain  deconvolution  experiments  in  sections  4.4.2  (two  parameter  method)  and  4.6.3.2 
(one  parameter  method).   Consequently,  the  time  domain  results  can  be  compared  to  the  two  frequency 
domain  methods,  and  vice-versa. 

The  input  and  output  waveforms  x(k)  and  y(k),  respectively,  are  shown  in  figure  5.14;  these  are 
the  signals  without  additional  noise  added  to  them.   Also,  the  time  scales  in  figure  5.14  are  expanded 
(50  ns./div.)  as  compared  to  those  in  figures  4.9  and  4.10  (200  ns/div.). 

Figure  5.15  shows  the  estimated  impulse  response  d(k)  using  the  filter  parameters  L  =  2,  S  =  5, 
A  =  3.15  applied  to  x(k)  and  y(k).   Compare  this  result  with  those  of  figures  4.7b  and  4.12  (2-parameter 
F.  D.  method)  and  figure  4.34  (1-parameter  F.D.  method). 

Figure  5.16  shows  the  estimated  impulse  response  d(k)  using  the  same  filter  L=2,  S=5,  A=3.15, 
but  with  noise  added  to  x(k)  and  y(k)  to  yield  SNR-k  and  SNR-y  values  of  40  dB  each.   Compare  this 
result  with  figures  4.7f  and  4.23  (2-parameter  F.D.  method)  and  figure  4.35  (1-parameter  F.D.  method). 
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Figure  5.1   The  time  positions  and  durations  of  x(k) ,  h(k) ,  y(k)  and  h^ (k) . 
x(k)  and  h(k)  are  two  given  sequences,  each  commencing  at  k  =  0; 
their  durations  are  N  and  M,  respectively.   y(k)  is  calculated 
by  convolving  x(k)  and  h(k),  and  the  result  commences  at  k  =  1. 
h^ (k)  is  h(k)  calculated  from  y(k)  and  x(k)  using  deconvolution. 
The  calculations  employ  eqs.  (5.2.1-8)  and  (5.2.2-6)  or  (5.2.2-21) 
for  convolution  and  deconvolution,  respectively. 


96 


CO 
I>Q- 


o 
> 

c 
o 
o 
<u 


1-1     f^ 

1-1 

C 


<U 

c 


o 


oo 

en 

\IO 

C 

g 

cn 

^la 

1-1 

1  CO 

(U 

iiJS 

3 
p 

*-« 

•H 

x: 
H 

X 

CO 

rH 

d 

00 


<u 

cfl 
•H 
O 

c 

<u 

x; 

4-1 

B 
o 
v^ 

<4-l 

o 
•u 

0) 

e 


d 

60 
•H 


97 


o 
> 

c 
o 
u 


re 

e 

60 


CO 


O 


U3  ^ 

W  ^— ' 

rt  ^>^ 

rH 

c 

(1)  to 

4-1 

<U  /^ 

O  ^— ' 

to  ^X 

-a  en 


c 

•H 
O 


ex: 

c 
o 
ei- 
m 

(U 

VI 

(U 

tn 
& 

e 


x; 


to 

c 


3 
Cu 

3 
O 

-O 

c 
to 


3 
C 


J-l 

D. 
3 
V- 
1-1 
O 

o 


o 

4-1 

e 


3 
00 

1-1 
O 


)-i 

to 

e 
o 
u 


to 

c 

00 
•H 
CO 

x: 

4-1 

o 

ja 

)-i 
o 


0) 
M 
3 
00 
■H 


98 


o 
> 

o 
u 

0) 

> 

-ri 
U 

u 

0) 


o 
o 


a: 


T3        TJ 
C 

01      rt 


>. 

^ 


o 


OJ 

to 
c 
o 
a. 
to 

01 

<i) 

to 


3 


0) 

3 


c 

•H 

tn 

OJ 
QJ 
iJ 

1 

a; 
to 

•H 

o 

c 

OJ 


O 

J-i 
OJ 

e 


Z>(f> 


o 
> 

c 
o 
o 

-o 

0) 
> 


(U 


CO 

o 


03 

c 
o 

■r- 
4J 
C3 

0) 

4-1 


QJ  O 

(-1  o 

o  o 

CO  r-l 


-a     -a 
c 

0)         CO 


c 


o 


c 


0) 

I 


oo 

0) 

C 

o 

o 

c 

Qi: 

QJ 

^UJ 

OJ 

£ 

1  1- 

Vj 

AJ 

en 
a. 

E 
o 
u 

CN 

in 

mcd 

E 

TJ 

0) 

o^^ 

•H 

O 

x; 

)-4 

-J 

»^ 

QJ 

■u 

00 

100 


c 
o 

•H 
4-> 
3 
iH 
O 
> 

c 
o 
o 

0) 

-a 


03 
C 


pa 

T3 


o 


CO     a: 
i-j     ^^ 

c 

0)        CO 
CJ        ^-^ 

en   ?X 


■o 


en 


0) 


6 
o 
u 

TD 
O 

JZ 
4-1 
0) 

e 


I 

LO 

OJ 
M 

00 
•H 


0) 
J-i 
3 
oo 


OO) 

4-1 

CO 

U 

"'g 

o 

to 

•H 

• 

ISIi-i 

•H 

cn 

cn 

CtKL 

4-1 

3 

o 

•H 

oo: 

s: 

Cu 

4J 

xm 

4J 

-1 

CO 

H- 

>^ 

O 

(U 

v-« 

^ 

4-) 

t3 

•H 

T3 

K 

O 

a> 

CO 

O 

to 

•H 
CO 

4J 

3 

4-1 

a 

J^ 

c 

• 

o 

T-l 

r-l 

m 

CO 

Q 

y— S 

t: 

c 

^« 

^ 

OJ 
4-1 

•H 

i-iO 

ex: 

a, 

m 

\ 

0) 

(n 

3 

"a 

^1- 

o 

o 

o 

CO 

^. 

en 

QJ 

^1 

2 

4-1 

o 

101 


01 
N 

•H 

u 


a; 


o 


M 

C 

o 

cn 

0) 


3 


2:yj<r 


en 

OJ 


U 

3 

00 


no 
O 

x: 
4J 

01 

e 

c 
o 


o 
> 
c 
o 
<J 


102 


T3 
(U 
N 

•H 
U 
CO 

tH 
3 
00 
0) 

u 


JO 
•a 

c 


O 


CO 

c 
o 
a. 
en 
tu 
^1 

q; 
en 

iH 

3 
a 


4) 

3 
GO 


O 
4J 

0) 

ca 
o< 

e 
o 
o 


o 

4-1 
0) 

e 

c 
o 


'^  o 

^  > 

13  O 
O 

0)  0) 

J-1  -a 

CO 

e  oj 

•U  ^ 

tn  u 

0)  nj 


3 
00 
•H 


103 


I 
I 


(a)      L-   Filter.      L  =   3, 

f        =   340  MHz.      Vertical 

CO 

scale  10  dB/Div.  Horizontal 
scale  50  MHz/Div. 


CO 


Pi 


(b)   S,  A-Filter, 


=  5,  A  =  1 
Vertical 


f   =  250  MHz. 

CO 

scale  10  dB/Div.  Horizontal 
scale  50  MHz. 


CO 


(c)  Combination  of  L  and  S 
Filters.  The  L  filter 
f   =  340  coincides  with 

CO 

the  S,  A  Filter  first  maximum. 
Vertical  and  horizontal  scales 
are  10  dB/Div.  and  50  MHz,  Div. 


Figure  5.9   Transfer  properties  L  and  S,A  filters  and  their  cascade  properties, 
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A  =  0.50 


A  =  1.00 


A  =  1.50 


A  =  1.80 


A  =  2.00 


A  =  2.05 


A  =  2.08 


A  =  2.40 


A  =  5.00 


Figure  5.10   Change  in  the  transfer  properties  of  the  S,  A 
Filter  as  A  is  increased.   S  =  5. 
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Figure  5.11   Examples  of  composite  filters. 
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Figure  5.12   Additional  examples  of  composite  filtering  scales  are 
the  sam.e  as  in  figure  5.12. 
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Figure  5.14   Input  and  output  waveforms  x(k)  and  y(k)  used  in  the  frequency 
and  time  domain  deconvolution  experiments  before  noise  is  added 
to  the  signals.   These  waveforms  are  identical  to  those  in 
figures  4.9  and  4.10  except  that  the  abscissae  are  50  ns/Div. 
(as  compared  to  200  ns/Div.). 
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Figure  5.15   Result  of  the  time  domain  deconvolution  d(k)  the  L2  and  S=5 ,  A=3.15 
filters  applied  to  the  same  waveforms  used  in  the  frequency  domain 
deconvolution  experiments  (figure  5.14).   Compare  to  the  frequency 
domain  deconvolution  result,  figure  4.12. 
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Figure  5.16   Result  of  time  domain  deconvolution  d(k)  when  the  input  and  output 
waveforms  are  corrupted  by  the  addition  of  35  and  6  dB  of  noise, 
respectively,  (figures  4.20  and  4.21)  to  yield  SNR-x  40  dB  and 
SNR-y  40  dB.   The  filters  applied  to  each  waveform  are  the  same  as 
those  in  figure  5.15.   Compare  this  result  to  the  frequency  domain 
result,  figure  4.23. 
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6.   SUMJIARY  AND  CONCLUSIONS 
6 . 1   Sununary 

The  objective  of  the  studies  reported  here  was  to  investigate  methods  of  deconvolution  to  be  applied 
to  time  domain  sequences  such  as  those  that  represent  the  signals  in  sampled  or  discrete  data  systems. 

In  the  preceeding  chapters  several  deconvolution  methods  have  been  explored  analytically  and  experi- 
mentally.  After  some  introductory  comments  and  some  general  presentations  of  the  deconvolution  problem 
and  error  criteria,  two  distinct  domains  for  deconvolution  of  time  domain  sequences  were  investigated: 
(1)  the  frequency  transform  sequence  domain,  or  frequency  domain,  and  (2)  the  time  sequence  domain,  or 
time  domain. 

Two  methods  of  discrete  frequency  domain  deconvolution  were  studied;  the  first  used  a  two  parameter 
filter  whose  parameters  were  subjectively  selected  by  the  operator  to  filter  the  transfer  function  H(n) 
prior  to  its  inversion  to  the  time  domain  to  yield  the  estimated  impulse  response,  d(k).   The  second 
method  was  one-parameter  optimal  technique  in  which  the  frequency  domain  filter  was  optimally  adjusted  to 
minimize  an  objective  function  consisting  of  the  weighted  sum  of  an  error-energy  and  roughness-energy 
(departure  from  smoothness).   The  two  methods  were  applied  to  the  measurement  of  the  impulse  response  of 
a  300  meter  long  coaxial  cable  for  which  experimental  input  and  output  time  domain  sequences  were  known; 
various  degrees  of  computer  generated  noise  (pseudo-random  sequences)  were  added  to  the  relatively  noise 
free  experimental  data  sequences  in  order  to  simulate  the  deconvolution  of  noisy  data.   Also,  in  addition 
to  the  coaxial  transmission  line  experiments,  the  two-parameter  method  was  applied  to  the  measurement  of 
the  impulse  response  of  a  broadband  antenna,  while  the  one-parameter  method  was  employed  to  obtain  the 
impulse  response  of  a  waveshaplng  filter. 

In  the  time  domain  studies,  two  methods  were  also  investigated,  discrete  classical  deconvolution 
and  discrete  iterative  deconvolution;  however,  the  latter  method  was  not  applied  to  any  experimental  data. 
The  classical  method  directly  inverted  the  discrete  convolution  equation  to  obtain  the  unknown  impulse 
response.   In  the  iterative  method,  the  discrete  convolution  equation  is  used  to  form  successive  approxi- 
mations to  the  impulse  response;  the  resultant  approximate  impulse  responses  are  convolved  with  the 
input  signal  to  yield  an  output  sequence  which  is  then  compared  to  the  known  output  sequence.   When  the 
two  sequences  agree  within  some  specified  limits,  the  iteration  process  is  stopped  and  the  N-th  iterated 
impulse  response  is  declared  to  be  equal  to  the  sought-after  or  unknown  impulse  response.   Both  methods 
were  applied  in  a  simulation  study  using  specified  x(k),  h(k),  and  y(k)  sequences  uniquely  related  to 
each  other  through  the  convolution  equation,  in  the  absence  and  presence  of  simulated  noise.   Also,  a 
simple  regularization  filter  was  used  to  demonstrate  the  improvement  provided  by  filtering  to  offset  the 
effects  of  noise.   The  classical  method  and  a  more  sophisticated  (composite)  filter  were  applied  to  the 
same  coaxial  cable  experimental  data  which  were  used  in  the  frequency  domain  deconvolution  studies; 
thus  the  results  of  the  frequency  and  time  domain  methods  could  be  compared  to  each  other.   The 
synthesis  and  properties  of  the  composite  discrete  time  domain  were  discussed  and  developed. 

6.2   Conclusions 

For  the  deconvolution  methods  considered  here  the  frequency  domain  one-parameter  method  provided 
the  smoothest  results  without  any  obvious  aberations  such  as  Gibbs-ef f ect.  The  two-parameter  method 
may  be  useful  in  some  interactive  filtering  situations,  but  for  the  most  part,  it  is  only  of  academic 
interest.  The  one-parameter  method  is  an  optimal  filtering  method,  and  as  such  the  weighting  between 
error  and  smoothness  (energies)  can  be  controlled. 

The  time  domain  deconvolution  methods  are  still  in  their  infancy  and  are  purely  experimental. 
Consequently,  much  work  remains  to  be  done  to  bring  them  to  fruition  so  that  they  can  be  easily  applied. 
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Nothing  encountered  in  this  investigation  rules  out  the  promise  held  for  time  domain  methods,  i.e. 
the  elimination  of  transform  filtering. 
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piled from  the  world's  literature  and  critically  evaluated. 
Developed  under  a  worldwide  program  coordinated  by  NBS  under 
the  authority  of  the  National  Standard  Data  Act  (Public  Law 
90-396). 


NOTE:  The  principal  publication  outlet  for  the  foregoing  data  is 
the  Journal  of  Physical  and  Chemical  Reference  Data  (JPCRD) 
published  quarterly  for  NBS  by  the  American  Chemical  Society 
(ACS)  and  the  American  Institute  of  Physics  (AlP).  Subscriptions, 
reprints,  and  supplements  available  from  ACS,  1 155  Sixteenth  St., 
NW,  Washington,  DC  20056. 

Building  Science  Series — Disseminates  technical  information 
developed  at  the  Bureau  on  building  materials,  components, 
systems,  and  whole  structures.  The  series  presents  research  results, 
test  methods,  and  performance  criteria  related  to  the  structural  and 
environmental  functions  and  the  durability  and  safety  charac- 
teristics of  building  elements  and  systems. 

Technical  Notes — Studies  or  reports  which  are  complete  in  them- 
selves but  restrictive  in  their  treatment  of  a  subject.  Analogous  to 
monographs  but  not  so  comprehensive  in  scope  or  definitive  in 
treatment  of  the  subject  area.  Often  serve  as  a  vehicle  for  final 
reports  of  work  performed  at  NBS  under  the  sponsorship  of  other 
government  agencies. 

Voluntary  Product  Standards — Developed  under  procedures 
published  by  the  Department  of  Commerce  in  Part  10,  Title  15,  of 
the  Code  of  Federal  Regulations.  The  standards  establish 
nationally  recognized  requirements  for  products,  and  provide  all 
concerned  interests  with  a  basis  for  common  understanding  of  the 
characteristics  of  the  products.  NBS  administers  this  program  as  a 
supplement  to  the  activities  of  the  private  sector  standardizing 
organizations. 

Consumer  Information  Series — Practical  information,  based  on 
NBS  research  and  experience,  covering  areas  of  interest  to  the  con- 
sumer. Easily  understandable  language  and  illustrations  provide 
useful  background  knowledge  for  shopping  in  today's  tech- 
nological marketplace. 

Order  the  above  NBS  publications  from:  Superintendent  of  Docu- 
ments, Government  Printing  Office,   Washington,  DC  20402. 
Order  the  following  NBS  publications — FIPS  and  NBSIR's — from 
the  National  Technical  Information  Services,  Springfield,  VA  22161 . 

Federal    Information    Processing    Standards    Publications    (FIPS 

PUB) — Publications  in  this  series  collectively  constitute  the 
Federal  information  Processing  Standards  Register.  The  Register 
serves  as  the  official  source  of  information  in  the  Federal  Govern- 
ment regarding  standards  issued  by  NBS  pursuant  to  the  Federal 
Property  and  Administrative  Services  Act  of  1949  as  amended. 
Public  Law  89-306  (79  Stat.  1127),  and  as  implemented  by  Ex- 
ecutive Order  11717(38  FR  12315,  dated  May  II,  1973)  and  Part  6 
of  Title  15  CFR  (Code  of  Federal  Regulations). 

NBS  Interagency  Reports  (NBSIR) — A  special  series  of  interim  or 
final  reports  on  work  performed  by  NBS  for  outside  sponsors 
(both  government  and  non-government).  In  general,  initial  dis- 
tribution is  handled  by  the  sponsor;  public  distribution  is  by  the 
National  Technical  Information  Services,  Springfield,  VA  22161, 
in  paper  copy  or  microfiche  form. 
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